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1.   INTRODUCTION 

One  of  the  most  important  areas  of  study  in  engineering  and  physics 
is  that  branch  of  statistical  mechanics  known  as  transport  theory.   A 
wide  variety  of  physical  phenomena  involve  particle  transport  processes. 
The  radiant  energy  transfer  in  a  stellar  atmosphere,  the  number  of 
particles  emerging  from  a  radiation  shield,  and  the  power  distribution 
in  a  nuclear  reactor  are  all  governed  by  the  distribution  of  particles 
which  interact  stochastically  with  the  atoms  of  the  host  medium.   The 
determination  of  such  distributions  is  the  central  problem  of  transport 
theory. 

The  solution  of  this  problem  becomes  much  more  difficult  when  the 
scattering  processes  involved  are  anisotropic.   Anisotropic  scattering 

is  encountered  in  a  wide  range  of  transport  phenomena  such  as  Compton 

4 
scattering  of  gamma  photons,   light  transmission  through  clouds  and 

hazes,  *  multigroup  neutron  transport  with  fine-energy-group 

7  8 

structures,   and  elastic  electron-nuclear  scattering. 

Scattering  anisotropy  can  result  from  several  factors.   Consider, 

for  example,  the  elastic  scattering  of  neutrons.   When  viewed  in  the 

center-of-mass  coordinate  system,  such  scattering  is  isotropic  for 

low-energy  neutrons.   At  neutron  energies  above  *\>0.1  MeV,  scattering 

generally  becomes  anisotropic,  with  the  scattering  cross  section  usually 

9 
being  forward  peaked. 

Scattering  anisotropy  Is  also  induced  by  the  transition  from  the 

center-of-mass  coordinate  system  to  the  laboratory  coordinate  system, 

which  is  the  frame  of  reference  most  often  used  in  transport 

calculations.   Scattering  which  is  isotropic  in  the  center-of-mass 


system  becomes  forward  peaked  in  the  lab  system.   While  this  effect  is 
not  significant  for  heavy  nuclei  (in  which  case  there  is  little 
difference  between  the  two  systems) ,  it  is  very  important  for  light 
nuclei. 

Use  of  the  multigroup  approximation  of  the  energy  variable  in  the 
transport  equation  introduces  group-to-group  scattering  (or  transfer) 
cross  sections  which  are  often  nonzero  over  only  a  limited  range  of 
scattering  angles.   This  limited  angular  support  is  a  direct  consequence 
of  the  constraints  for  energy  and  momentum  conservation.   Such 
anisotropy  becomes  more  severe  as  finer  multigroup  structures  are  used. 

The  traditional  method  of  treating  anisotropic  scattering  in 
transport  calculations  has  been  to  expand  the  angular  dependence  of  the 
scattering  cross  section  in  a  finite  series  of  Legendre  polynomials. 

Such  an  expansion  is  mathematically  convenient  to  use  because  the 

2 
Legendre  polynomials  obey  an  "addition  theorem   which  allows 

significant  analytical  simplification  of  the  scattering  source  term  in 

the  transport  equation.   This  method  is  accurate  if  a  sufficiently  high 

order  expansion  is  used.   However,  a  highly  anisotropic  cross  section 

often  requires  a  very  high  order  expansion  to  represent  the 

cross  section  accurately.   Use  of  a  prematurely  truncated  Legendre 

expansion  often  introduces  spurious  oscillations  in  the  approximated 

cross  section.   These  oscillations  in  turn  can  produce  unrealistic 

fluctuations  in  the  calculated  angular  fluxes,  and  may,  in  fact,  even 

lead  to  estimates  of  negative  angular  fluxes. 

The  failure  of  the  Legendre  expansion  method  to  produce  physically 

realistic  angular  fluxes  occurs  only  in  problems  which  are  characterized 


by  highly  anisotropic  angular  fluxes  as  well  as  by  highly  anisotropic 
scattering.   Although  the  degree  of  anisotropy  in  the  angular  flux  is 
seldom  known  a  priori,    problems  with  anisotropic  sources  and  vacuum 
boundary  conditions  as  well  as  anisotropic  scattering  can  be  expected  to 

yield  anisotropic  flux  distributions,  especially  near  the  sources  and 

23 
boundaries . 

To  avoid  the  problem  of  negative  fluxes,  many  alternative  methods 

of  cross  section  approximation  have  been  proposed.   One  of  the  simplest 

remedies,  which  is  known  as  the  transport  approximation,  is  to  replace 

an  anisotropic  scattering  cross  section  that  is  highly  peaked  in  the 

forward  direction  by  an  isotropic  component  and  a  forward-scattered 

component.    The  forward-scattered  component  is  most  simply  represented 

by  a  delta  function.   Although  this  transport  approximation  can  provide 

good  accuracy  for  problems  which  are  characterized  by  only  a  small 

12 
degree  of  anisotropy,   it  often  gives  poor  results  for  highly 

11,12 
anisotropic  problems. 

13 
Razani   proposed  a  modified  transport  approximation  in  which 

singly-scattered  particles  are  accounted  for  exactly  and  the  effect  of 

higher-order  scattering  is  treated  by  the  transport  approximation.   He 

found  that,  for  radiation  transport  through  homogeneous  slabs,  the  error 

in  the  transport  approximation  was  considerably  reduced  by  using  this 

12 
modified  transport  approximation.   Lathrop   examined  an  extended 

transport  approximation,  obtained  by  adjusting  the  coefficients  of  a  P 

(i.e.  ,  first  order  Legendre)  approximation.   He  applied  this  method  to 

several  neutron  transport  problems  and  found  it  to  be  more  accurate  than 


the  delta-function  transport  approximation,  but  not  as  accurate  as 

14 
higher  order  standard  Legendre  expansions.   Bell  et  at.        examined  the 

extended  transport  approximation  for  higher  order  Legendre  expansions 

and  concluded  that  it  was  an  effective  method  for  neutron  transport 

problems. 

Attia  and  Harms   used  a  partial-range  Legendre  polynomial 
expansion  of  the  scattering  cross  section.   While  this  representation 
yields  more  accurate  cross  section  fits  than  a  standard  full-range 
expansion,  it  does  not  lend  itself  readily  to  use  in  discrete-ordinates 
transport  codes.   Pearlstein   proposed  an  expansion  of  the  scattering 
cross  section  in  terms  containing  quadratic  Bessel  functions.   He  found 
that  such  an  expansion  accurately  modeled  scattering  cross  sections  for 
a  variety  of  elements  with  fewer  terms  than  a  standard  Legendre 
expansion  required.   However,  the  use  of  quadratic  Bessel  functions  is 
mathematically  inconvenient,  and  it  has  not  been  established  whether 
such  an  expansion  would  be  feasible  in  discrete-ordinates  calculations. 

Carter  and  Forest   utilized  a  step  function  representation  of  the 

scattering  cross  section  in  Monte  Carlo  transport  calculations. 

18 
Takeuchi   assumed  a  step  function  approximation  of  the  scattering  cross 

section  with  respect  to  the  lethargy  mesh,  and  a  step  function 

approximation  of  the  angular  flux  with  respect  to  both  the  angular  mesh 

and  the  lethargy  mesh.   These  approximations  were  utilized  in  the  PALLAS 

computer  code,  which  numerically  solves  the  integral  form  of  the 

transport  equation. 

A  more  direct  approach  to  avoid  the  use  of  Legendre  polynomial 

cross  section  expansions  is  to  use  a  transfer  matrix  whose  elements  are 


the  exact  cross  sections  for  transfer  from  one  discrete  direction  to 

19 
another.   While  this  possibility  had  been  considered  in  the  1960's,   it 

was  not  fully  developed  at  that  time  due  to  the  associated  requirement 

20 
of  large  computer  memories.   Odom,   who  finally  implemented  this 

technique  in  plane-geometry,  discrete-ordinates  neutron  transport 

calculations,  referred  to  it  as  the  "exact  kernel"  method.  He  found 

that  the  exact  kernel  method  provides  accurate  angular  fluxes  even  for 

highly  anisotropic  problems. 

21 
Mikols   reduced  the  computational  effort  of  calculating  exact 

kernel  cross  sections  by  assuming  a  triangular  representation  of 

group-to-group  neutron  transfer  cross  sections.   He  also  developed  the 

"order  of  angular  coverage"  concept  for  determination  of  the  minimum 

order  of  numerical  quadrature  required  for  accurate  transport 

calculations. 

22 
flyman   applied  the  exact  kernel  technique  to  discrete-ordinates 

gamma  photon  transport  problems.   He  showed  that  the  exact  kernel 

technique  yields  far  more  accurate  results  for  such  problems  than  does 

the  standard  Legendre  expansion.   He  also  developed  a  semi -analytical 

technique  for  rapid  evaluation  of  exact  kernel  cross  sections. 

23 
Hong   applied  the  exact  kernel  method  to  neutron  inelastic 

scattering.   He  developed  an  exact  kernel,  discrete-ordinates  transport 

code  which  is  applicable  to  plane,  spherical,  cylindrical,  and  two-angle 

plane  geometries .   He  also  introduced  a  method  to  evaluate  the 

scattering  source  term  using  piecewise  polynomial  approximations  for  the 

angular  flux  and  the  transfer  cross  section. 


While  the  exact  kernel  technique  offers  greater  accuracy  than  the 
Legendre  expansion  method,  it  has  the  disadvantage  of  being  incompatible 
with  standard  discrete-ordinates  transport  codes.   Because  these 
standard  codes  almost  always  are  based  on  Legendre  polynomial 
expansions,  they  would  require  modification  in  the  scattering  source 
term  calculation  in  order  to  utilize  exact  kernel  cross  sections. 
Beranek  and  Conn   suggested  a  method  by  which  a  discrete  transfer  cross 
section  expansion  can  be  utilized  in  standard  discrete-ordinates 
transport  codes.   This  method  involves  the  generation  of  pseudo-Legendre 
expansion  coefficients  which  will  reproduce  the  exact  kernel  cross 
sections  for  scatter  between  any  discrete  directions  of  the  quadrature 
set  used.   However,  the  number  of  expansion  coefficients  required  to 
reproduce  all  the  exact  kernel  cross  sections  increases  rapidly  with  the 
quadrature  order,  rendering  this  method  impractical  for  highly 
anisotropic  problems. 

In  this  work,  the  method  suggested  by  Hong  for  the  semi-analytical 
evaluation  of  the  scattering  source  term  is  developed  more  fully  and  is 
applied  to  several  transport  problems.   This  new  method  is  more  accurate 
for  highly  anisotropic  problems  than  is  the  conventional  Legendre 
expansion  method.   In  addition,  the  new  method  eliminates  the  problem  of 
angular  coverage  encountered  in  exact  kernel  transport  calculations. 

In  Chapter  2  a  general  development  of  the  multigroup,  discrete- 
ordinates  transport  equations  is  presented.   Evaluation  of  the 
scattering  source  terra  is  discussed  for  both  the  Legendre  expansion 
method  and  the  exact  kernel  method.   In  Chapter  3,  the  approximated 
scattering  kernel  method  is  developed.   This  semi-analytical  evaluation 


of  the  scattering  source  term  utilizes  piecewise  polynomial  expansions 
of  the  transfer  cross  section  and  the  angular  flux.   Use  of  these 
polynomial  expansions  allows  the  integration  of  the  scattering  source 
term  to  be  performed  analytically,  thereby  obviating  the  need  for 
numerical  quadrature. 

In  Chapter  4,  the  approximated  scattering  kernel  method  is  applied 
to  several  slab  albedo  problems.   The  results  are  compared  to  those 
obtained  using  the  exact  kernel  and  Legendre  expansion  methods. 
Finally,  the  conclusions  reached  from  these  problems,  as  well  as 
suggestions  for  further  study,  are  presented  in  Chapter  5. 


2.   SOLUTION  OF  THE  TRANSPORT  EQUATION  BY  THE  DISCRETE-ORDINATES  METHOD 

2.1  The  Transport  Equation 

The  transport  of  neutral  particles  (e.g.,  neutrons  and  photons) 
through  a  host  medium  is  described  by  the  linearized  Boltzmann  equation. 
This  equation,  which  neglects  particle-particle  interactions,  is  a 
linearized  form  of  the  equation  derived  by  Boltzmann  more  than  a  century 
ago  in  his  study  of  the  kinetic  theory  of  gases.   The  derivation  of  the 
linearized  Boltzmann  equation  (hereafter  referred  to  as  "the  transport 
equation")  can  be  found  in  many  texts  (see,  for  example,  Refs.  1  and  2) 
and  will  not  be  repeated  here. 

In  this  work,  we  are  concerned  only  with  steady-state  transport 
through  a  homogeneous,  non-multiplying  slab  for  cases  in  which  the 
particle  flux  density  possesses  azimuthal  symmetry.   Thus  we  seek 
solutions  of  the  steady-state,  one-dimensional  transport  equation,  which 
can  be  written  as 

II  -g-;  *(x,E,U,t)  +  a(E)*(x,E,u,i|i)  =  Q(x,E,u,t|i) 

■     1      2tt 

dE'    du'    dili'  a  (I,+E,8*'JS)  *(x,E'  ,u'  ,  +  ' )  ,     (2.1) 
0    >-l    h 


where 


u  =  the  cosine  of  the  polar  angle  between  a  particle's  velocity 
vector  and  the  positive  x-axis  (see  Fig.  2.1), 

x  =  the  distance  travelled  into  the  slab, 

E  =  the  particle's  energy, 

i>   =  the  azimuthal  angle  between  the  z-axis  and  the  projection  of 
the  particle's  velocity  vector  onto  the  slab  face, 


Scatter 
Part  id 


Fig.  2.1.   The  coordinate  system  for  particle  transport  in  a  slab. 

The  slab  dimensions  in  the  Y  and  Z  directions  are  assumed 
to  be  infinite. 
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$  =  the  angular  flux  density, 
a  =  the  total  macroscopic  cross  section, 
Q  =  extraneous  (flux-independent)  sources, 
a  =  the  macroscopic  scattering  cross  section, 
and 

Q  =  a  unit  vector  in  the  direction  of  particle  travel. 

The  medium  is  assumed  to  be  isotropic  so  that  a  is  rotationally 
invariant  (i.e.,  a  depends  on  w  =  ft?"fi,  not  on  Q1  and  U   individually). 

In  terms  of  the  incident  particle  direction  H1  (v'»<|if)  and  final 
particle  direction  R(u»$)i  the  cosine  of  the  scattering  angle  is  given 

by2 

u  =  Sf«£  =  uyf  +  (l-p2)*2  (l-vi'2)^  cos*  (2.2) 

where  $  is  the  difference  in  the  azimuthal  angles  of  the  vectors  ft'  and 
S  (i.e.,  <f>  =  i|»f-ijj). 

For  azimuthally  symmetric  problems,  one  often  solves  for  the 
azimuthally  integrated  flux  density.   Integration  of  Eq.  (2.1)  over  the 
azimuthal  angle  i^  yields  the  transport  equation  for  the  azimuthally 
integrated  flux  density,  viz.: 

U   3^  *(x,E,u)    +  a(E)*(x,E,u)    =  Q(x,E,p) 

(1  (-2TT 

dE'        dy'  diji  Os(E'*E,ui)    *(x,E*,u'),  (2.3) 


where 


f2it 

■Kx.E.u)    =         d<l>  *(x,E,y,i|i)  (2.4) 

0 
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and 

Q(x,E,u)  =    I   d*  Q(x,E,u,i|0  .  (2.5) 

0 

Due  to  the  dependence  of  oj  on  cos<J>  =  cos(iJj '-<|0  ,  the  azimuthal  integral 
over  all  f  on  the  right  hand  side  of  Eq.  (2.3)  can  equivalently  be 
replaced  by  an  integral  over  all  i/i T  or  $.   Because  we  are  concerned  only 
with  the  change  in  azimuthal  angle  upon  scattering,  we  will  find  it  most 
convenient  to  work  with  the  variable  <j> ,  and  thus  Eq.  (2.3)  is  rewritten 


u  t—  *(x»E,u)  +  ff(E)*(x,E,y)  =  Q(x,E,u) 


C00      r  1      j-  2lT 

+   dE'   dp*    d$  a    (E'+E.w)  *(x,E',u')  •  (2.6) 

0-10 

Even  in  this  simplified  form,  the  transport  equation  is  much  too 
complex  to  solve  analytically  for  realistic  geometries.   Thus  it  is 
common  practice  to  introduce  further  approximations  such  as  discretizing 
the  energy,  angular,  and/or  spatial  variables.   By  using  such 
approximations,  one  is  left  with  sets  of  algebraic  equations  rather  than 
an  integrodifferential  equation.   The  resulting  algebraic  equations  lend 
themselves  readily  to  numerical  solution  techniques. 

2.2  Energy  Approximations  for  the  Transport  Equation 
2.2.1  The  Multigroup  Transport  Equations 

The  energy  dependence  of  Eq.  (2.6)  is  commonly  approximated  by 
employing  the  multigroup  method.   This  formulation  is  accomplished  by 
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discretizing  the  energy  variable  into  G  contiguous  energy  groups,  as 
shown  schematically  below: 

Group  g 

I 1 1 v— H 1 1 \ — l 1 1 

EG    EG-1   EG-2       Eg+1  Eg    Eg-1       E2    El    E0 


ENERGY 

The  standard  multigroup  convention,  which  will  be  followed  in  this  work, 

is  to  index  the  energy  groups  such  that  the  group  index  g  increases  with 

decreasing  energy.   The  energy  level  E  lies  within  group  g  if 

E   <  E  <  E 
g     -  g-1 

To  derive  the  multigroup  equations,  we  first  define  the  group 
angular  flux  density  and  group  extraneous  source  as 

*  (x,u)  =  dE  4>(x,E,u)  (2.7) 

E 
g 

and 

fEg-l 
Q  (x,u)  =  dE  Q(x,E,u)  .  (2.8) 

E 
g 

Upon  integrating  Eq.  (2.6)  over  energy  group  g  and  replacing  the 
integration  over  energy  in  the  scattering  source  term  by  a  sum  of 
integrals  over  all  the  energy  groups,  we  obtain  the  following  multigroup 
equations: 
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u  •£-;  «  (x,u)  +  0  *  (x,u)  ■  Q  (x,y) 

G 

+   I  S  ,   (x,u)   ,    g  =  1 G  (2.9) 

g'=l  8  g 

where  the  group  total  cross  section  is  defined  by 

Vi 

dE  o(E)  *(x,E,u) 
E 

a  E  — §— ; ,  (2.10) 

8      -Vl 


dE  *(x,E,u) 

E 
8 


the  scattering  source  is  defined  by 


1     2ir 


S  .   (x,u)  E   du'  I   d*  a  .   (u>)  *  '(x,u')  ,  (2.11) 

g  ""-g        J     I  g  *%  g 

-1    0 

and  the  group-to-group  scattering  cross  section  is  defined  by 


E„  i    A , 


g-1 
dE 

E 


g'-l 

dE'  a    (E'+E.oi)  *(x,E',u) 


E 


a  ^  (a>)  i  — S 1 .         (2.12) 

1 
dE'  4(x,E*,u) 


'  g'-l 


E 
g 


The  evaluation  of  the  group  cross  sections  through  the  use  of  Eqs. 
(2.10)  and  (2.12)  presupposes  knowledge  of  the  angular  flux  density, 
which  is  not  known  until  the  transport  equation  itself  is  solved.   This 
problem  is  usually  circumvented  by  assuming  that  the  energy  dependence 
of  the  flux  is  separable  from  the  spatial  and  angular  dependence,  i.e., 
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*(x,E,u)  ■  W(E)  f(x,p).   With  this  assumption,  Eq.  (2.10)  can  be 
rewritten  as 


where 


1  8-1 

dE  0(E)    W(E) 

>E 
g 

0       — 

g 

A 
g 

A       5 

g 

rVl 

dE  W(E)    . 

>E 

g 

(2.12) 


(2.13) 


In  a  similar  manner,  Eq.  (2.12)  can  be  rewritten  as 

1     r   g-1   (■  g'-l 
.(u)  =  j^t  dE       dE'  W(E')  o-g(E'+E,w)  .        (2.14) 


g  +g 


E       E  , 
g       g 


For  very  fine  energy  meshes,  the  weight  function  W(E)  is  often 
assumed  to  be  constant  over  each  energy  group.   For  broader  energy 
meshes,  a  variety  of  approximations  have  been  used.   For  example,  in 
neutron  transport  in  fission  reactors,  one  often  uses  a  fission  spectrum 

for  W(E)  in  the  MeV  range,  a  1/E  spectrum  in  the  epithermal  range,  and  a 

25 
Maxwellian  spectrum  in  the  thermal  energy  region. 

2.2.2  The  One-Speed  Transport  Equation 

In  many  particle  transport  problems,  the  particles  are  assumed  to 
be  characterized  by  a  single  energy.   This  assumption  is  a 
simplification  of  the  multigroup  equations  in  that  only  one  energy  group 
need  be  utilized.   The  resulting  equation  is  known  as  the  "one-speed"  or 
"one-group"  transport  equation.   It  is  also  referred  to  as  the  "constant 
cross-section  approximation",  since  if  the  cross  sections  are  postulated 


15 


to  be  independent  of  energy,  the  one-speed  equation  results.   Although 
this  treatment  may  seem  rather  crude,  it  actually  has  a  great  deal  of 
practical  application,  *   and  is  the  basis  of  most  analytical  studies 
of  the  transport  equation. 

The  one-speed  transport  equation  can  be  derived  quite  simply  from 
the  multigroup  equations.  Let  us  first  write  the  multigroup  equations 
in  the  following  form  (cf.  Eqs.  (2.9)  and  (2.11)): 


M  -r-  *  (x,u)  +  a  4  (x,u)  =  Q  (x,y) 

»x  s       s  s       s 


G   »     2* 


g 


I      [  du*  I  d*  a  ,  M    *  .(x.u*)  ,  g  -  1,...,  G.   (2.15) 
i  =i  j      J  8   8      g 


Since  we  are  considering  only  one  energy  group,  we  can  drop  the  group 
subscripts  and  delete  the  summation  over  the  incident  energy  groups  in 
the  scattering  source  term  to  obtain  the  one-speed  transport  equation, 
viz.  : 

a  f1    l2v 

n|-«(x,u)  +  CTt(x.u)  ■  Q(x,u)  +   du'    d$  <JaU)    4>(x,y').  (2.16) 

-1     0 

It  is  often  convenient  to  recast  the  transport  equation  into 
dimensionless  form.   To  do  so  we  express  distances  in  terms  of  the 
collision  mean-free-path  length,  also  referred  to  as  the  "optical 
thickness".   We  define  this  dimensionless  distance  as 


(2.17) 


so  that 


\h 


\-=o\-.  (2.181 

3x     3z 


Eq.  (2.16)  can  now  be  rewritten  as 

uer  y-*(s,w)  ♦  CF*(*,U)  -  OQ(l.u)  ♦  I  dw'    d*  as(u)  *(z,u'),     (2.19) 


where 

*(z,u)  =  *(x(z),u)  (2.20) 


and 


Q(x(z)  ,u)  (2.21) 


Q(z.y)  = 


Equation  (2.19)  can  be  simplified  by  expressing  the  differential 

scattering  cross  section  a  (to)  as  the  product  of  the  total  scattering 

cross  section  a  and  a  scattering  distribution  function  f  (<i>) ,  viz.: 
s  "  s 

(J  (a)  S  0  f  (id)  ,  (2.22) 

s       s  s 

where  f„(n))  is  normalized  to  unity,  i.e., 

dQ  f  (tt)  >  1  .  (2.23) 

s 

4ir 

Substitution  of  Eq.  (2.22)  into  Eq.  (2.19)  yields 

3  f1     f  2,t 

uo  |-  *(z,n)  +  o*(z,u)  =  oQ(z,u)  +  ag   du  '    d*  f  (u)  »(z,u')   (2.24) 

-1     0 
or,  in  dimensionless  form, 
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f1  f2lT 

u  T~  «(z,u)    +  4(z,u)    ■  Q(z,u)    +  c       du'  d<t>   fs(o))    4(z,u')    ,  (2.25) 

-1  0 

where  c,  the  "mean  number  of  secondary  particles  emitted  per  collision," 
is  given  by 


c  =  —  .  (2.26) 


2.3  The  Discrete-Ordinates  Equations 
2.3.1  Angular  Discretization 

The  angular  dependence  on  the  transport  equation  is  commonly 
approximated  by  the  discrete-ordinates  method.   In  this  method,  the 
angular  variable  u  is  discretized  into  a  set  of  N  discrete  directions 
{u.}  at  which  the  angular  flux  is  to  be  evaluated.   The  scattering 
source  term  is  evaluated  by  numerical  quadrature  where  the  \i .  values  are 
the  quadrature  ordinates.   The  set  of  corresponding  quadrature  weights 
is  denoted  by  {w.}.   With  this  approximation,  Eq.  (2.9)  can  be  rewritten 


"i  k  *g(x'pi)  +ag  VX*V  =  Qg(x>Mi) 


G 
+     J     S    ,      (x,u.)    ,    i  =   1,...,   N;      g  =    I,...,   G  (2.27) 


where 


V^'V  E  J,  wj     d*  V*g(B)  V^'V  •  (2>28) 


N  f2TT 

I 

J-l 
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Equations  (2.27)  and  (2.28)  represent  one  form  of  the  discrete-ordinates 
transport  equations . 

The  evaluation  of  the  group-to-group  scattering  source  S  f   (x,u.) 
of  Eq.  (2.28)  is  dependent  upon  the  method  used  to  express  the 
group-to-group  scattering  cross  section.   Two  of  these  methods,  the 
Legendre  polynomial  expansion  method  and  the  exact  kernel  method,  will 
be  examined  in  sections  2.4  and  2.5.   A  new  method  of  evaluating  the 
scattering  source,  which  is  the  subject  of  this  work,  will  be  examined 
in  Chapters  3  and  4.   This  alternative  method  is  based  on  direct 
evaluation  of  Eq.  (2.11)  by  analytical  integration. 

The  accuracy  which  can  be  obtained  in  solving  the  discrete- 
ordinates  equations  when  Eq.  (2.28)  is  used  to  evaluate  the  scattering 

2 
source  term  is  largely  dependent  on  the  quadrature  set  used.    In 

general,  one  would  like  to  use  a  set  which  is  large  enough  to  describe 

adequately  the  angular  detail  in  the  fluxes,  yet  small  enough  that 

excessive  computational  effort  is  not  required.   The  choice  of  such 

an  optimum  set  is  typically  problem  dependent,  especially  when 

anisotropic  scattering  is  involved.   Failure  to  choose  an  appropriate 

quadrature  set  can  lead  to  serious  errors  in  the  calculation  of  the 

angular  fluxes. 

Unfortunately,  there  is  no  standard  procedure  for  choosing  a  priori 
an  adequate  set.   The  choice  is  usually  made  on  the  basis  of  experience 
or  trial  and  error.   However,  as  a  general  rule,  the  following  criteria 
should  be  met: 

(1)  Projection  Invariance.      For  one-dimensional  slab  geometry  with 
azimuthal  symmetry,  the  discrete  directions  {u.}  should  be  symmetric 
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about  u  =  0.   However,  if  one  knows  that  the  angular  flux  is  peaked  near 
a  certain  direction  u, ,  it  may  be  advantageous  to  tailor  a  nonsymmetric 
quadrature  set  with  several  discrete  directions  clustered  near  u.. 

(2)  Positivity  of  the  scalar  flux.     The  scalar  flux 

r1  N 

•  (x)  =    du  *  (x,u)  =  I     w  *  (x,u.) 

g        J  6  j_i    1   S      x 

-1  X     l 

should  be  always  positive.   Choosing  the  w  >  0  will  ensure  positivity 
of  the  scalar  flux  (provided  the  angular  flux  values  are  positive) . 

(3)  Accurate  Evaluation  of  Angular  Integrals.     The  scalar  flux  and 
scattering  source  should  be  evaluated  accurately  with  a  minimum  of 
quadrature  ordinates  and  weights. 

Two  commonly  used  quadrature  sets  for  one-dimensional  geometries 
are  the  Gaussian  quadrature  set  and  the  Lobatto  quadrature  set.   Values 
of  the  ordinates  and  weights  of  these  quadrature  sets  for  various  values 

of  N  can  be  found  in  Ref.  29.   Because  the  angular  flux  is  discontinuous 

2 
at  a  plane  interface  at  u  ■  0,   odd  order  quadrature  sets  are  not  used 

in  plane  geometry  transport  calculations  (since  they  contain  a 

quadrature  ordinate  at  y  =  0) .   In  order  to  avoid  this  possible 

discontinuity,  one  can  split  up  the  angular  integration  range  into  two 

parts,  -1  <  u  <  0  and  0  <  u  <  1,  and  perform  Gaussian  quadrature 

separately  over  each  sub-range.   This  approach  is  known  as  the  double  P 

(or  DP„)  method.   In  plane-geometry  transport  problems,  the  DP   method 

has  often  been  found  to  give  numerical  results  which  are  superior  to 

2 
those  obtained  from  a  standard  full-range  expansion. 
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2.3.2  Spatial  Discretization 

For  some  simple  cases  (e.g.,  a  one-group  problem  with  isotropic 
scattering  and  no  extraneous  sources) ,  the  discrete-ordinates  equations 
(2.27)  can  be  solved  analytically.   This  approach  was  used  by  Wick  and 

Chandrasekhar  in  the  original  development  of  the  discrete-ordinates 

25 
method.    However,  for  most  realistic  transport  problems,  Eqs.  (2.27) 

are  too  difficult  to  solve  analytically  and  hence  must  be  solved 

numerically. 

One  method  of  effecting  a  numerical  solution  of  Eqs.  (2.27)  is  to 

form  a  set  of  finite-difference  equations  by  discretizing  the  spatial 

variable  x  into  a  set  of  spatial  nodes  {x,  }.   Let  the  left  boundary  be 

denoted  by  x1  and  the  right  boundary  by  x^+,  (see  Fig.  2.2).   The 

spatial  derivative  of  the  flux  is  then  approximated  by  a  finite 

difference  scheme,  viz.: 


|-»  K+u  V  s    s  a 8  ■  (2-29) 

3x  g  k+4,  i  A^ 


where 


and 


*k+h  '  2 


\+h   "  *k+l    *k  • 

We  thus  obtain  the  finite  difference  form  of  the  multigroup 
discrete-ordinates  equations  as 


(2.30) 


(2.31) 


a^  ^Vv^i'  -  *g(xk'u±)]  +  °g  Vvv^i' 
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VWl'  '   k 


.,  K;   i  -  1, 


N;   g  =  1,...,  G 


(2.32) 


where  the  group  total  source  q  is  defined  by 

G 

VWi*  '=  VwV  +  Jml  WwV 


(2.33) 


Fig.  2.2.   Spatial  discretization  of  a  one-dimensional  slab. 


2.3.3  Solution  of  the  Spatially  Discretized  Discrete-Ordinates 
Equations 

Before  the  set  of  equations  represented  by  Eq.  (2.32)  can  be 

solved,  it  is  necessary  to  reduce  the  number  of  unknowns  by  assuming  a 

relation  between  the  cell-edged  fluxes  and  cell-centered  fluxes.   The 

most  commonly  used  method  is  the  diamond  difference  scheme  which  uses 
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i^Wi 


M.J  0 


(2.34) 


Substitution  of  Eq.  (2.34)  into  Eq.  (2.32)  yields 


t  (a^ .Wj),  k-  X K;   i-1 N;  g  =  1,...,  G 


(2.35) 


which  represents  K  x  N  x  G  equations  in  K  x  N  x  G  unknowns,  *  (x,  ,U.). 
(The  incident  flux  densities  at  the  outer  surfaces  of  the  slab,  x.  and 


x_  .  ,  are  assumed  known  from  the  boundary  conditions.) 

Equation  (2.35)  can  be  solved  for  *  (x,  ,»".)  in  terms  of 

yvv  as 


ywv 


i  _  gA+is 


lv. 


k+'-i 


i  ♦-4-S4 

2U± 


(2.36) 


or  for  *  (x,  ,U.)  in  terms  of  *  (x,  ,,11.)  as 

8   K   1  g   K+l   1 


i  »  tVi 


2U, 


,(xk'V 


•g^i-V-  -rf  y=wV 


g  k+'5 
2u. 


(2.37) 


These  two  results  permit  the  evaluation  of  the  angular  flux  densities  at 
all  the  spatial  nodes  by  starting  at  one  of  the  slab  surfaces  (where  the 
incident  flux  is  known)  and  sweeping  inward  through  the  spatial  mesh 
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along  the  direction  of  particle  travel.   Equation  (2.36)  must  be  used 
for  u.  >  0  and  Eq.  (2.37)  for  u   <  0.   This  procedure  minimizes  the 
accumulation  of  roundoff  errors. 

The  number  of  spatial  nodes  and  the  number  of  discrete  ordinates 
are  not  independent  of  each  other.   From  examination  of  Eqs.  (2.36)  and 
(2.37),  it  is  seen  that  to  ensure  the  positivity  of  the  left-hand  side 
of  each  equation,  the  following  condition  must  be  met: 


S  k+% 


2*i 


<  1  (2.38) 


Thus,  the  maximum  cell  width  A,  ,  is  constrained  by  the  smallest  cosine 

of  the  polar  angle  (i.e.,  the  value  of  u.  closest  to  zero).   As  more 

discrete  directions  are  used,  the  minimum  |u.  |  becomes  smaller,  and 

hence  the  number  of  spatial  nodes  required  for  positivity  of  the 

computed  fluxes  increases. 

The  solution  of  the  discrete  ordinates  equations  is  complicated  by 

the  fact  that  the  source  term  q  (x,  x  »u.)  is  dependent  on  the  unknown 

flux  $    (x.  ,  ,u  )  (except  for  a  purely  absorbing  medium,  in  which  case 
g  k+*5  i 

a    ,   (ui)  is  zero)  .   The  usual  method  of  solving  the  discrete  ordinates 
g  ''•g 

equations  is  thus  to  employ  an  iterative  scheme  in  which  the  source  term 
is  approximated  more  accurately  with  each  iteration.   The  procedure  is 
as  follows: 

(1)  Estimate  the  initial  source  term  q  (x.  .  ,u  )  at  each  spatial 
cell  midpoint  using  any  reasonable  source  distribution. 

(2)  Solve  for  the  angular  flux  densities  using  Eqs.  (2.36)  and 
(2.37)  as  appropriate. 
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(3)  Calculate  the  cell-centered  fluxes  using  Eq.  (2.34). 

(4)  Reconstruct  the  source  terms  using  Eqs.  (2.33)  and  (2.28). 

(5)  Repeat  steps  (2)  -  (4)  until  the  computed  flux  densities 
converge. 

The  rate  of  convergence  of  this  iterative  process  depends  on  both 
the  nature  of  the  scattering  involved  and  the  thickness  of  the  slab. 
For  the  extreme  case  of  no  scattering  (a  =0)  the  method  converges  in  a 
single  iteration.   As  the  amount  of  scattering  increases,  so  in  general 
will  the  number  of  iterations  required  for  convergence.   For  optically 
thick  slabs  with  values  of  the  scattering  ratio  c  close  to  unity, 

i    i    30 
convergence  is  extremely  slow. 

Because  of  the  need  to  reduce  computational  time  (and  hence  cost) 
in  discrete  ordinates  calculations,  various  methods  have  been  developed 
to  accelerate  the  convergence  of  the  iterative  routine.   Some  of  the 
more  common  methods  are  reviewed  in  Ref.  31.   The  method  used  in  this 
work  is  a  simple  overrelaxation  (SOR)  scheme  which  is  implemented  as 
follows . 

Consider  the  calculation  of  the  scattering  source  terms 
S  ,   (x.  .i  u.)  for  iteration  v  +  1.   These  values,  denoted  by 

S  ,   ft  ..u.),  can  be  expressed  as  a  vector  functional  of  the 

20 
previous  iterate  scattering  source  values,   i.e. 


s^g  Wi' =  '(te'Vwv}]  •  (2-39) 
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20 
It  has  been  found   chat  convergence  may  be  accelerated  by  calculating  a 

modified  source  term,  viz. : 


-  ^{ww^}  (2-4o) 

where  [I]  is  the  identity  matrix  and  [6  ]  is  a  diagonal  matrix  defined 
by 

[6V]  =  dv[I].  (2.41) 

The  quantity  d  is  a  scalar  between  0  and  1  which  may  vary  with  the 
iteration  number. 

Ryman  obtained  good  acceleration  results  by  setting  d  initially 
to  some  small  value  (=0.1)  and  multiplying  it  by  a  constant  (sl.l)  with 
each  iteration  until  an  upper  limit  (s0.2)  was  reached.  This  procedure 
was  found  to  provide  good  results  in  this  work  as  well. 

Having  developed  the  multigroup  discrete-ordinates  equations  and 
described  the  method  by  which  they  are  solved,  we  now  turn  our  attention 
to  the  evaluation  of  the  scattering  source  term. 

2.4  Evaluation  of  the  Scattering  Source  Term  by  the  Legendre  Expansion 
Method 

The  accuracy  which  can  be  obtained  in  the  solution  of  the  discrete 

ordinates  equations  is  dependent  on,  among  other  factors,  the  accuracy 

obtained  in  the  evaluation  of  the  scattering  source  term.   In  Section 

2.3.1  it  was  mentioned  that  the  evaluation  of  the  scattering  source  term 

is  dependent  upon  the  method  used  to  express  the  group-to-group 
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scattering  cross  section.   In  the  next  two  sections,  we  will  examine  two 
of  these  methods. 

The  most  widely  used  method  is  to  represent  the  angular  dependence 
of  the  scattering  cross  section  by  a  truncated  Legendre  polynomial 
expansion,  i.e. , 


'••«■«  "  £  (tt)  °*.g'*8  V»>  ■  <2-^ 


where  the  expansion  coefficients  are  given  by 

1 

"n.g'+g  =  2*  I  du  Wu)  V"0  (2,43) 

-1 

and  P.  is  the  Legendre  polynomial  of  order  %. 

25 
Application  of  the  addition  theorem  for  Legendre  polynomials   to 


Eq.  (2.42)  yields 


'  J0  W\   °*.8'*8  (P<(1°  P*(W,) 


g  -"g 


£ 


2i  Tiiyfp>)p>,)cos(m*)}'  (2-44) 

m=l  ' 


where  the  P.  are  associated  Legendre  polynomials  of  the  first  kind. 
Upon  insertion  of  Eq.  (2.44)  into  Eq.  (2.28),  the  terms  containing 
cos(m<f>)  vanish  upon  integration  over  0.   The  result  is 


WX'V  -T  ,lx  wjJ0(2£+"  "*.■•*«  '^  W  V(s"V-  (2-45> 
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Equation.  (2.45)  can  also  be  expressed  in  terms  of  the  angular  flux 
moments  as 


where  the  moments  of  the  flux  are  given  by 

,(x)  =  2lt  f  du  »  ,(x,u)  P.(u).  (2.47) 


Jt.S 

-1 

While  Eqs.  (2.45)  and  (2.46)  provide  convenient  methods  for  the 
evaluation  of  the  scattering  source  term,  they  suffer  from  a  serious 
drawback.   A  low  order  expansion  of  two  to  nine  terms  (typically  used  in 
reactor  physics  calculations)  may  fail  to  model  the  scattering  cross 
section  adequately  .   This  is  particularly  true  when  highly  anisotropic 
scattering  is  involved.   Such  anisotropic  scattering  is  common  in 
neutron  transport  problems  involving  elastic  scattering  from  light 
elements  and  with  very  fine  energy  group  structures.   As  an  example, 
consider  the  hydrogen  multigroup  scattering  cross  section  of  Fig.  2.3. 
The  eighth-order  Legendre  expansion  is  seen  to  exhibit  oscillatory 
behavior,  with  the  cross  section  expansion  actually  having  negative 

values  over  portions  of  the  w  range.   The  use  of  such  an  expansion  can 

20 
lead  to  the  calculation  of  oscillatory  and  even  negative  flux  values. 

While  the  use  of  higher  order  expansions  could  mitigate  this  problem, 

the  extra  computational  effort  involved  often  makes  such  an  approach 

unfeasible. 

In  spite  of  this  drawback,  the  Legendre  expansion  method  is  widely 

used  and,  in  most  situations,  provides  satisfactory  results.   For  many 
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situations  in  which  the  scattering  cross  sections  would  appear  to 
require  a  very  accurate  angular  description,  low-order  Legendre 
expansions  are  often  justified  if  the  angular  flux  is  not  highly 
anisotropic  or  if  only  angularly  integrated  quantities  (e.g.,  scalar 
fluxes)  are  of  interest.   Only  for  those  problems  in  which  both  the 
angular  flux  and  the  scattering  cross  sections  are  highly  anisotropic  is 
it  necessary  to  resort  to  a  more  accurate  method,  such  as  the  exact 
kernel  technique. 

2.5  Evaluation  of  the  Scattering  Source  Term  by  the  Exact  Kernel  Method 

For  problems  in  which  the  Legendre  expansion  method  produces 
inaccurate  or  nonphysical  fluxes,  more  accurate  results  can  be  obtained 
by  using  a  non-expansion  method.   This  method,  referred  to  by  Odom  as 
the  "exact  kernel"  method,  involves  the  use  of  a  transfer  matrix 
composed  of  scattering  cross  sections  for  every  u   to  u .  transfer.   Use 
of  the  exact  kernel  method  guarantees  positivity  of  the  computed  fluxes 

(provided  Eq.  (2.38)  is  satisfied)  and  has  been  shown  to  provide 

20  21,23  22 

accurate  results  for  both  neutron        and  gamma  photon   transport 

problems. 

The  exact  kernel  form  of  the  scattering  source  term  is  simply 

where  a  ,   (u.-*u.)  is  the  exact  kernel  cross  section  for  transfer  from 
g  *S     J   i- 

polar  direction  u  to  polar  direction  u . .   Comparison  of  Eqs.  (2.28)  and 
(2.48)  shows  that  the  exact  kernel  cross  section  is  defined  as 
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tj.  (M.-y.)  =   d*  a  ,   (a)(y,,U.,<fO)  ,  (2.49) 

*g  J   i    J     g  *g    1  J 


where  u(|i,u',$)  is  given  by  Eq.  (2.7).   Due  to  the  dependence  of 

a  ,   (u>)  on  cos((i,  which  is  symmetric  about  tt  in  the  interval  [0,  2it]  , 

Eq.  (2.49)  can  also  be  written  as 


a  .   (u.-nO  =2   d*  a  ,   (ui).  (2.50) 

g '*g  j  i      J     g  *g 

0 

Several  techniques  have  been  developed  to  facilitate  the  evaluation 
of  Eq.  (2.50).   These  techniques  include  restriction  of  the  azimuthal 
integration  range  and  the  use  of  low-order  piecewise  polynomial 
approximations  for  a    ,   (id)  .   The  method  used  to  evaluate  exact  kernel 
cross  sections  for  the  transport  calculations  performed  in  this  work  is 
presented  in  Sec.  3.2.1.   For  more  information  on  the  evaluation  and  use 
of  exact  kernel  cross  sections,  the  interested  reader  is  referred  to 
Refs.  22  and  23. 

Although  the  exact  kernel  method  does  provide  more  accurate  results 
in  problems  with  a  high  degree  of  anisotropy  than  does  the  Legendre 
expansion  method,  it  does  have  drawbacks.   The  most  obvious  of  these  is 
the  problem  of  cross  section  storage.   In  the  Legendre  expansion  method, 
the  only  cross  section  terms  required  are  the  L  +  1  expansion 
coefficients  for  an  L'th  order  cross  section  expansion.   In  the  exact 
kernel  method,  however,  it  is  necessary  to  store  an  N  x  N  matrix  of 
cross  section  values  (N  being  the  number  of  discrete  ordinates) .   For 
transport  problems  involving  a  large  number  of  discrete  ordinates  and/or 
many  energy  groups,  it  is  easy  to  see  that  the  number  of  exact  kernel 
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cross  sections  required  becomes  very  large.   While  this  storage 
requirement  is  a  disadvantage  of  the  exact  kernel  method,  it  has  become 
less  of  a  problem  as  larger  computers  and  fast  access  peripheral  storage 
devices  have  been  developed. 

Another  problem  associated  with  the  exact  kernel  method  is  that  of 
angular  redistribution  of  scattered  particles.   This  problem  can  best  be 
illustrated  by  considering  a  highly  anisotropic  scattering  cross  section 
which  is  nonzero  only  over  a  small  portion  of  the  scattering  range 
a)  e  [-1,1]  (cf.  Fig.  2.3).   For  such  a  cross  section,  many  of  the  exact 
kernel  cross  sections  a(u.*u,)  will  be  identically  zero.   Whether  or  not 
a  particular  exact  kernel  cross  section  is  nonzero  depends  upon  the 
spacing  between  adjacent  u,  values  of  a  discrete-ordinates  quadrature 
set.   If  the  spacing  between  these  values  is  too  large,  particles 
traveling  in  direction  u.  will  never  scatter  into  other  directions. 
This  inability  of  particles  to  redistribute  angularly  remains  even  after 
multiple  scatters . 

In  order  to  ensure  that  angular  redistribution  can  occur,  it  is 

necessary  to  have  at  least  one  nonzero  value  of  cr(u  -*-u  )  for  some  i  4   j. 

27 
For  this  condition  to  be  met,  Mikols   has  shown  that  there  must  exist 

at  least  one  polar  quadrature  ordinate  8  =  cos   u.  such  that 

S  ,   <  6   <  B    .  (2.51) 

min  —  i  —  max 

The  values  of  6  ,   and  6    are  dependent  upon  the  nonzero  range  of  the 
min      max 

scattering  cross  section,  as  well  as  upon  the  initial  particle  direction 

6.  3  cos   u..   In  particular,  consider  a  cross  section  which  is  nonzero 

J         J 

only  over  the  ranee  [oi  .  ,  cu   ].   The  angles  8  .   and  8    are  given  by 
J  ■    mm   max         °     min      max  * 


min 


maxr(6.  -  6   ),  01  ,    8.  >  0  .   , 

j    max  *      j  —  mm 

6  .   -  8.  ,    8.  <  0  .   . 

min    j  j    min 
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(2.52) 


and 


(2.53) 


min [(6  +0   ),tt]   ,   6.  +  8  .   <ir, 
a      J       J    max  j    min  — 

max   I  2TT  -  e  -  e  .      ,e.+eJ>Tr, 

I  j    min  j    min 

where  0    3  cos   to  .   and  0  .   =  cos   u 

max         min      rain         max 

The  fundamental  concern  when  choosing  a  discrete-ordinates 
quadrature  set  for  exact  kernel  transport  calculations  is  whether  or  not 

the  set  can  "completely  model"  scattering  transfers  within  a  group  and 

27 
to  the  next  lower  group.    The  term  "completely  model"  refers  to  the 

ability  of  a  quadrature  set  to  produce  at  least  one  nonzero  value  of 

o(u  +u  ),  with  i  ^  j,  for  all  u.  values  of  the  quadrature  set.   This  is 

referred  to  as  first-order  coverage  of  a  quadrature  set.   In  the  same 

manner,  n'th-order  coverage  is  that  coverage  which  permits  transfer  from 

u.  to  at  least  n  different  u.  values.    In  general,  a  quadrature  set 

which  provides  first-order  coverage  is  sufficient  for  transport 

calculations.   While  larger  quadrature  sets  (with  a  higher  order  of 

coverage)  may  provide  some  increase  in  accuracy,  they  do  so  at  the  cost 

of  increased  computational  effort.   On  the  other  hand,  if  the  quadrature 

set  is  so  coarse  that  not  even  first-order  coverage  is  achieved,  the 

exact  kernel  method  yields  very  poor  results. 

In  the  next  chapter,  a  new  method  of  evaluating  the  scattering 

source  term  is  introduced.   This  new  method  avoids  the  angular  coverage 

problem  of  the  exact  kernel  method,  and  also  provides  better  accuracy 

for  highly  anisotropic  problems  than  does  the  Legendre  expansion  method. 
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3.   DEVELOPMENT  OF  AN  APPROXIMATED  SCATTERING  KERNEL  FOR 

SEMI-ANALYTICAL  EVALUATION  OF  THE  SCATTERING  SOURCE  TERM 

In  the  previous  chapter,  two  methods  for  evaluating  the  scattering 
source  term  of  Eq.  (2.28)  were  summarized  —  the  Legendre  expansion 
method  and  the  exact  kernel  method.   For  problems  characterized  by 
highly  anisotropic  fluxes  and  scattering  cross  sections,  the  Legendre 
expansion  method  can  lead  to  the  calculation  of  inaccurate  and  even 
unphysical  fluxes.   The  exact  kernel  technique  provides  much  more 
accurate  results,  but  can  suffer  from  angular  redistribution  problems. 

In  the  next  two  chapters  we  will  examine  a  method  based  on  the 
analytical  evaluation  of  the  scattering  source  term  of  Eq.  (2.11).   This 
method,  which  we  will  refer  to  as  the  "approximated  scattering  kernel" 
method,  will  be  seen  to  provide  better  accuracy  than  the  Legendre 
expansion  method,  and  to  alleviate  the  angular  redistribution  problem  of 
the  exact  kernel  method  for  highly  anisotropic  problems. 

3.1  Approximations  for  the  Scattering  Source 

In  Chapter  2  it  was  shown  that  the  group-to-group  scattering  source 
for  a  direction  u  is  given  by 

1      tZV 

0 


S  ,.  (y)  -   dp'   d*  a  ,.  (co)  4  (u-),      (3.1) 
g  *g      1  J     g  +g     g 


where  the  spatial  variable  has  been  dropped  for  simplicity  of  notation. 
Due  to  the  dependence  of  0  ,   (id)  on  cos4,  which  is  symmetric  about  1   in 
the  interval  [0,  2»] ,  Eq.  (3.1)  can  be  rewritten  as 

1     it 

S  ,.  (u)  =  2  [  dU'    d4  a    (w)  *  ,(U').         (3.2) 
g*g        I  1  g  *g     g 

-1    0 
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In  order  to  evaluate  Eq.  (3.2)  analytically,  it  is  first  necessary  to 

assume  a  functional  form  for  both  the  group-to-group  scattering  cross 

section  and  the  angular  flux  density.   In  this  section,  we  will  examine 

approximations  which  often  represent  the  shape  of  a  .   (w)  and  $  ,(u') 

g  +g        g 

more  accurately  than  the  traditional  Legendre  polynomial  expansions. 

3.1.1  Approximation  of  the  Scattering  Cross  Section 

In  previous  studies  of  transport  calculations  involving  anisotropic 
scattering,  it  has  been  shown  that  the  angular  range  of  the  scattering 
cross  section  a   ,      (id)  is  naturally  divided  into  subregions  bounded  by 

at  most  four  breakpoints  in. ,  i  =  1,  . . . ,  4,  such  that 

21-23 
-1  <  U,  <   u,  <  <d,  £  <o,  <  1.       The  values  u1  and  to,  represent  the 

minimum  and  maximum  cosine  of  the  scattering  angle,  respectively,  for 

which  scattering  is  permissible  within  the  constraints  of  energy  and 

momentum  conservation  for  transfer  from  group  g'  into  group  g.   Thus  the 

scattering  cross  section  a  .   (10)  is  zero  for  CD<td,  or  co>to, .   For 
g  -*-g  14 

anisotropic  scattering,  the  interval  [oo,  ,  id.  1  over  which  a    ,      (oi)  is 

1   4'  g'+g" 

nonzero  is  often  only  a  small  portion  of  the  entire  range  [-1,1]  (cf. 
Fig.  2.3). 

In  general,  the  breakpoints  in  the  scattering  cross  section  can  be 

23 
expressed  as 


*j  =  [[-1,  S(Eg,  Eg,_1),  1]]  ,  (3.3) 


m2   =  [[-1,  S(Eg_1,  Eg,_1),  1]]  ,  (3.4) 


[[-1,  S(Eg,  Egl),  1]]  ,  (3.5) 


and 
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u4  =  [[-1,  S(E   ,  E  ,),  1]]  ,  (3.6) 

where  the  notation  [[a,x,b]]  =  max[a,min(x,b) ] ,  i.e.,  =  x  if 
a<x<b,  =  b  if  x  >  b,  or  =  a  if  x  <  a.   The  function  S(E,E'),  which 
is  obtained  from  scattering  kinematics,  gives  the  relationship  among  the 
initial  and  final  particle  energies  and  the  cosine  of  the  scattering 

angle.   For  elastic  scattering  of  neutrons  from  a  nucleus  of  mass  number 

2 
A,  this  relationship  is  given  by 


0)  =  S(E.E')  E  i  (A+l)  f\r   -  (A-l)  /|H  .       (3.7) 


25 
while  for  Cotnpton  scattering  of   photons 

2  2 

m  c  ra  c 

a  =  S(E.E')    E    1   -  -2—  +  -%j-  ,  (3.8) 


2 
where  m  c  is  the  rest  mass  energy  of  the  electron  (0.511  MeV) . 

Under  certain  conditions,  the  four  breakpoints  defined  by  Eqs. 

(3.3)  -  (3.6)  are  not  all  distinct.   For  energy  group  structures  equally 

spaced  in  lethargy  [u  =  iln(E   /E  )]  or  Compton  wavelength 

2 
[X  =  (m  c  /E)  ] ,  to„  and  oj_  coalesce.   For  such  a  group  structure  we  can 

define  three  distinct  breakpoints  as 


Dj  =  [[-1,  S(E  ,  Eg, .p.  HI  .  (3.9) 


and 


,)2  =  [[-1,  S(Eg_1,  Eg,_1),  1]]  =  [[-1,  S(E  ,  E  ,),  11]  ,     (3.10) 


u3  =  [[-1,  S(E    ,  E  ,),  1]1  .  (3.11) 
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For  within-group  scattering  (e.g.,  one-speed  transport  problems),  ui9, 

(ii_,  and  ui.  are  coincident  at  u  •  1,   We  can  thus  define  two  distinct 
3       4 

breakpoints  as 


and 


Oj  =  [[-1,  S(E  ,  E  , _1),  1]]  ,  (3.12) 


U2  -  1  .  (3.13) 


The  utility  of  the  in.  breakpoints  lies  in  the  fact  that  the  shape 
of  the  scattering  cross  section  is  usually  very  smooth  in  each  m 

subrange  and  can  generally  be  well  modeled  by  piecewise,  low-order 

23 
polynomials  between  the  breakpoints.    Such  piecewise  polynomial  fits 

are  generally  far  superior  to  traditional  full-range  Legendre  polynomial 

expansions  (see,  for  example,  Fig.  2.3).   In  particular,  piecewise 

linear  fits  have  been  shown  to  be  a  good  approximation  for  neutron 

elastic  scattering  cross  sections  involving  fine-energy-group  structure 

or  scattering  from  light  elements,  while  piecewise  quadratic  fits  have 

been  shown  to  provide  good  results  for  Compton  scattering  of  gamma 

22 
photons.    In  this  work  we  will  consider  both  piecewise  linear  and 

piecewise  quadratic  models  for  the  scattering  cross  section.   We  can 

represent  a  general  model  for  both  cases  as 

k 

max  . 


g  +8  k^1      g  *g 


where 


.                        a.  ui  +  b.  id  +  c.       ,      -l<oi,<iii<u),     ,<1 

k        .   .        JTc  k           k                 -  Tc  -      -     k+1  -       ,,    ... 

an'^M    =  \  (3.15) 

0  ,      otherwise 


g  ->g 
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and  k    is  one  less  Chan  the  number  of  distinct  breakpoints  in  the 
max 

scattering  cross  section. 

Given  the  values  of  the  cross  section  at  the  breakpoints,  it  is  a 
simple  matter  to  compute  the  coefficients  of  Eq.  (3.15)  for  a  linear 
fit.   In  the  case  of  a  quadratic  fit,  it  is  necessary  to  know  the  cross 
section  at  one  additional  point  in  each  subdomain.   A  convenient  point 

to  use  is  the  midpoint  of  each  subdomain. 

22 
For  a  quadratic  fit  the  coefficients  are  given  by 

a^.  =  D  /D  ,  (3.16) 

bk  =  D2/D  ,  (3.17) 


and 


ck  =  D3/D  ,  (3.18) 


where 


D  =  \[ttkth  -  \*1]   "  Wis  "  "k+l1  +  «W"k  "  VV  '   OA9) 

Dl    =    °k!uk+%   "   <Vn]    "    ak+!s[uk   "   <"k+l]    +    Vl!\    "   "W     '  O-20) 

D2  =  "k^k-Mj  "  °k+l]    _  "iMijt'k  "  ak+l]    +  "W°k  -  <W    '        (3-2n 
and 

D3    =    "k^+Jj  °k+l    "   "k+1    <W    -   "W"k   <Vl    -    <\+l    °k] 

+  vi[\  ak+!i  -  %*  ak]  •  (3-22) 

In  Eqs.  (3.19)  through  (3.22)  u.  ,  is  the  cosine  of  the  scattering  angle 

at  the  midpoint  of  the  k'th  subdomain  and  a,  *   a,    i  ,  and  a,  ,  are  the 

k   k+%       k+1 

cross  section  values  at  oj,  ,  co.  !  ,  and  uv  ,»  respectively. 

For  a  piecewise  linear  fit  the  coefficients  of  Eq.  (3.15)  are  given 
by 


■k 
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0  ,  (3.23) 


=  JSil fe  ,  (3.24) 

k    ID,   ,  -  ID, 

k+l    k 


and 


(3.25) 


k+l 


3.1.2  Approximation  of  the  Angular  Flux 

The  second  approximation  needed  to  evaluate  the  scattering  source 
term  of  Eq.  (3.2)  is  a  representation  of  the  angular  flux  density. 
Perhaps  the  simplest  approximation  is  to  represent  the  angular  flux  by  a 

series  of  straight  lines  between  the  discrete  ordinates  (v.).   This 

2 
approximation  was  the  basis  of  Carlson's  original  S  method.    In  this 

work,  we  will  utilize  this  approximation  and  also  generalize  it  to 

higher-order,  piecewise  polynomial  fits. 

In  particular,  consider  a  set  of  (NM  +  1)  discrete  ordinates  {u.}, 

where  y  =  -1  and  U   ,   =  1.   We  can  divide  this  set  into  M  segments  of 
1  NM+ 1  ■ 

N  +  1  discrete  ordinates  each.   Let  the  ordinates  at  the  boundaries  of 

each  segment  be  denoted  byu.,j=l,...,M+l,  where  M ,  =  U,.  ,%„  ,. 

J  :     (j-l)N+l 

In  each  of  the  M  segments,  the  angular  flux  can  be  approximated  by 
an  N'th  order  polynomial  through  the  N  +  1  ordinates  in  that  segment. 
In  particular,  the  flux  in  the  j'th  segment  can  be  approximated  by  an 
N'th  order  Lagrange  polynomial  as 


jN+1 
*  00  "    I  L  (u)  *  (u  )    ,   i,  <  p  <  »'  (3.26) 

8      i=jN-N+l  l     8  *         J        J  ' 

where 
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jN+1  u-u 

L,(u)  =  n    —       ,   u.  <  u  <  h    .  (3.27) 

1      k=jN-N+l  Vk        3  "   "  J  +  l 


The  angular  expansion  of  Eq.  (3.26)  may  be  extended  to  the  whole  range 
u  e  [-1,11  as 

NM+1 

*JV)   -   I  L.  (u)  *  (ii  )    ,   -1  <  U  <  1        (3.28) 
8       i-1  1     8  i  -   - 

where  the  full  range  Lagrange  polynomials  are  defined  by 
M  I   jN+1   u-u 

Li(u)  E  I        n     in*   H([vVrVryi!)  •       (3>29) 

j=l  k=jN-N+l  Wi  Uk      1  J   3 
I  Mi        J 

In  Eq.  (3.29)  H(x)  is  the  unit  step  function  [i.e.  H(x)  =  0,  x  <  0; 

H(x)  =  1,  x  >  0],  so  that  L.  (y)  vanishes  unless  y  and  y.  belong  to  the 

j  '  th  segment . 

3.2  Decomposition  of  the  Source  Integrals 

Having  developed  explicit  piecewise  polynomial  approximations  for 
the  scattering  cross  section  and  the  angular  flux  density,  we  now 
consider  the  resulting  form  of  the  scattering  source  term  of  Eq.  (3.2). 
With  the  piecewise  quadratic  approximation  for  o  ,  (ut)   of  Eqs.  (3.14) 
and  (3.15),  the  scattering  source  term  can  be  written  as 

f     (Tr    max  , 
Sg,+gfti)  -  2  J  du'  I  d*  I     a*   (u)  i    ,(»')  (3.30) 

-1    0   k=1 

or,  equivalently , 
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k 
max  , 

h'+M    "     I      S  (ij)  '  (3-31) 


g  *g 


k-1 


where 


Sk(u)  E  2  |  du*  |  d*  <Jgi  +  g(u)  *g.(u').  (3.32) 

-1    0 
The  integration  of  Eq.  (3.32)  can  be  performed  analytically  using 

the  piecewise  polynomial  approximations  for  a  ,   (aj)  and  $  ,  (u ' ) .   In 

order  to  carry  out  the  integration,  it  is  first  necessary  to  decompose 

the  integration  ranges  of  u'  and  $  into  subranges,  in  each  of  which  the 

integrands  assume  a  distinct  functional  form. 

In  Sec.  3.1.1  it  was  noted  that  a  ,   ((d)  is  zero  outside  the 

interval  [to,  ,  oj,  .  ] .   Because  dj  is  a  function  of  u»  v',  and  <|>  Tcf.  Eq. 

(2.2)],  there  exist  corresponding  limits  of  the  u'  and  4>  integrals  of 

Eq.  (3.32)  outside  of  which  the  integrand  is  identically  zero.   In  this 

section  we  will  consider  the  specific  limits  of  the  y'  and  1(1  integrals 

for  Eq.  (3.32)  and  the  form  of  the  integrands  in  each  of  the  resulting 

integration  subranges. 

3.2.1  Restriction  of  the  Azimuthal  Integration  Range 

In  order  to  evaluate  the  inner  integral  (i.e.,  the  integral  over  (j>) 
of  Eq.  (3.32),  it  is  first  necessary  to  determine  the  precise  limits  of 
$  for  which  a    ,   (oj)  is  nonzero.   The  desired  limits  can  be  found  from 
the  iso-(u  contours  in  the  (u',  i)>)  integration  plane  as  shown  in  Figs. 
3.1  -  3.3.   These  contours  are  obtained  by  plotting  f  as  a  function  of 
u'  for  fixed  values  of  u  and  oi,  using  Eq.  (2.2).   The  minimum  and 
maximum  values  of  <j>  for  which  a  ,   (a))  is  nonzero  depend  on  u  and  u'  and 
equal,  respectively,  the  $  values  on  the  ui,  .  and  u),  contours  (or  their 
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limiting  values  of  $  ■  0  and  <)>  =  it)  .   The  $  integration  range  is  thus 

k+1         k 
reduced  from  [0,  Til  to  [if    (y.U1),  <t>  (y,y')],  where,  from  Eq.  (2.2), 


k           -1         \  "  W 
Au.y')  =  cos  l[[-l,  £ 5-p  ,  1]1.         (3.33) 

(l-liZ)'5  (1-u'V 


The   inner  integral  of   Eq.    (3.32)    thus  becomes 

!k+l 


k, ,    =    I*  k        ,...,    _    |T      ,M    ,_      2 


I   (y,y')    =       d*   a    ,+    (u)    =  d<(>    (a^     +  b^  +  cfe)  ,        (3.34) 

0  <t>k" 

which,   upon  use  of  Eq.    (2.2),    can  be   integrated   to  give 

Ik(y,y')    =    [aku2u'2  +  bkMi'    +  cfc]    (*k  -  f>k+1) 

+   (l-y2)!s(l-u'2)I's[2alcuu'   +  bkl  (sini(.k  -  sin*k+1) 

+  Y~   (1-U2)(l-P'2)[*k  -  *k+1  +  sin*k  cos<t>k  -  sin*k+1  cos*k+1].    (3.35) 


Substitution  of  Eq.  (3.34)  into  Eq.  (3.32)  yields 
1 
Sk(u)  =  2  j  dy'  Ik(y,y')  4>  ,(y').  (3.36) 

-1 
Before  we  proceed  with  the  evaluation  of  the  y'  integral  of  Eq. 

(3.36),  it  is  instructive  to  examine  the  scattering  source  as  defined  by 

Eqs.  (3.31)  and  (3.36),  i.e., 

k      i 
max  A 

S<,,_(")  -  I  2   dy'  I*(y,y')  *  ,(y').         (3.37) 
8  8      k=1   I  8 

This  expression  could  be  evaluated  by  numerical  quadrature  to  yield 

max  NM+1     . 

V+e(u)  =   I       I  2  n»,ll  )  »  »,(,.),      (3.38) 

8   8        k-1  j=l  J    J   g    J 
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Comparison  of  Eq.  (3.38)  with  Eq.  (2.48)  shows  that  the  exact  kernel 
cross  section  defined  in  Section  2.5  can  thus  be  calculated  as 

k 
max  . 
Tk 


WW  =  2  Jx  T  Cui»V  (3-39) 


3.2.2  Evaluation  of  the  Polar  Integral 

In  order  to  evaluate  the  u'  integral  of  Eq.  (3.36)  analytically,  it 
is  first  necessary  to  decompose  the  u'  integration  range  into  several 
subranges,  in  each  of  which  I  (y,u')  and  $  ,(u')  assume  a  single 

functional  form.   The  form  of  I  (u,u')  in  each  subrange  is  dependent  on 

k      k+1 
(J)  and  0   *  which  take  on  different  forms  depending  on  u1  .   From  Figs. 

k      k+1 
3.1  -  3.3,  it  is  seen  that  $  and  <j>    assume  a  single  functional  form 

—  1  2  h  2 

(namely,  0,  tt,  or  cos   { (w-uu? )/[  (1-y  )  (1-y1  )]})  in  at  most  five 

subintervals  of  y'  e  [-1,1],   Let  the  u'  breakpoints  which  define  these 

subintervals  be  denoted  by  iC,,*  t  -   1,...,  6.   These  breakpoints,  which 

are  dependent  on  p,  ov»  and  w,  .,  are  the  y1  values  at  which  the  to.  and 

u,  .  contours  in  the  u'-<t>  space  of  Figs.  3.1  -  3.3  intersect  the  lines 

$  *  0  and  (J)  =  it.   In  particular,  u'   =  -1,  yj      =   1,  and  the  other  four 

values  are  the  ordered  roots  of  Eq.  (2.2)  with  0  =  0  or  0  =  tt,  i.e., 

{uk2'Uk3'Uk4'Uk5}  "  °rdered  Kl'Ks'K^'^  '  (3-40) 

where 

^2  =  M\  +  [(l-V2)(l-^)]h   ,  (3.41) 

U^3  "  W\  "  [(1-U2)(1-«^)]'S  ,  (3.42) 

»U~»\+i  +  [(^'"-Vi'''1  ■  "-43> 
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and 

^5  =  ^k+i-  raVm-^)]*  .  (3.44) 

k      k+1 
Unless  m.  or  to,  ,  is  equal  to  +u  ,  $   and  <J>    will  both  equal  0  or 

tt  in  the  intervals  [u*  ,  p'  ]  and  [u*  ,  u^ ft ]  -   Thus  I  (u,u!)  will  vanish 

in  these  intervals,  so  they  can  be  neglected  when  evaluating  Eq.  (3.36). 

Further  decomposition  of  the  u'  integration  range  can  be  achieved 

by  considering  the  piecewise  polynomial  representation  of  the  angular 

flux.   In  each  [u,jU,  .]  segment  the  flux  is  represented  by  Eq.  (3.26), 

where  the  interpolating  polynomial  of  Eq.  (3.27)  can  be  expressed  more 

simply  as 

N 
L.(u)  =  I     C^   un   ,   S  <  u  <  u  (3.45) 

n=0  J        J 

Because  the  coefficients  C  are  different  in  each  of  the  j  segments,  the 

integrand  of  Eq.  (3.36)  assumes  a  different  form  in  each  segment. 

There  are  thus  6  breakpoints  (v^p}  in  the  u'  integration  range  of 

Eq.  (3.36)  imposed  by  the  functional  form  of  I  (u,yT)»  and  M+l 

breakpoints  {u.}  imposed  by  the  flux  interpolating  polynomials. 

However,  u,'   ■  U.  =  -1  and  u'   =  uw  .  ■  1,  so  the  number  of  distinct 
Kl     1  ko    M+l 

breakpoints  is  at  most  M+5  (note  that  it  Is  possible  for  some  of  the 
other  u'   and  u.  breakpoints  to  be  coincident).   These  breakpoints 
divide  the  u'  integration  range  into  at  most  M+4  subranges,  in  each  of 

which  the  integrand  of  Eq.  (3.36)  assumes  a  single  functional  form. 

k  f- 

Let  tv.)  represent  the  ordered  values  of  (iC?^  anc*  tv.}.   Equation 

(3.36)  can  then  be  rewritten  as 
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k 
.        M+4  A+l 
Sk(y)  -  2  I  dy'  Ik(y,y')  *  ,(y').  (3.46) 

\ 

Use  of  Eqs.  (3.28)  and  (3.45)  allows  Eq.  (3.46)  to  be  expressed  as 


.         NM+1   M+4  N   ... 
Sk(y)  -  2  I       \l       I     FkUn(y)U   (y  ),  (3.47) 


1=1  u=i  n-o        8  1 


where 


k 

Fki±n(y)  SCM  '  dy'  y'n  Ik(y,y').  (3.48) 

n  J  k 

The  problem  of  evaluating  the  scattering  source  components  S  (y)  is 

thus  reduced  to  the  evaluation  of  the  integrals  F  (y)  as  defined  by 

Eq.  (3.48).   In  the  next  section  we  will  examine  the  analytical 
evaluation  of  these  integrals. 

3.3  Evaluation  of  the  F   n  Integrals 

In  order  to  evaluate  the  integrals  F    (y)  of  Eq.  (3.48),  we  must 

consider  the  exact  form  of  the  integrand  in  each  of  the  integration 

k   k 
subranges  [v  ,  v        } .      Substitution  of  Eq.  (3.35)  into  Eq.  (3.48)  yields 

FkJlin(y)  =  Ci{GUn(y;k)  -  Gkin(y;k+1)} ,  (3.49) 


where 

k 

„k«.n,  ,  ,,  -  (  +  ,  ,  /  ,n.   2,2   ,    ,     ,  ,k' 
G   (y;k')  =      dy'  jy'  [a^v   y'  +  b^y1  +  ck]  ♦ 


+  y'n(l-y2)!i(l-y'2)!i[2akyy'  +  bk]  sin  $k' 

+  yiy'n  (l-y2)(l-y'2)[<t.k'  +  sin^'cos*1"' ]}.   (3.50) 
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The  exact   form  of  G        (u;k')    depends   on  the   form  taken  by  if      .      In 

k' 
Section  3.2   it  was   shown   that   0        can,    in  general,    take   on  three  forms 


(cf.    Eq,    (3.33)),   although  only  one   in  a  given  integration  subrange. 

kin 
Thus  G        (u;kT)    is   given  by  one  of    three  possible   forms,  viz.: 

,k' 


Case    1: 


0; 
Gkin(]i;k')    =  0    . 


(3.511 


Case   2: 


xk' 


Jtta,     .,,  JfV2        \a-u2)}(,k     ,n+3        .k.n+3 

G      (y;k')  ■»    b I(n73T-    <Vi>        -  <V 


b,.U 


kM    I,  k      xn+2        .  k.n+2] 


n+2      (Vl' 


fTc  ak(1"y2)H,  k      ,n+l         ,  k.n+11' 


[n+l  2(n+l)    J  | 


M' 


(3.52) 


k'               -if       V    "  vv 
Case  3:        +       =  cos       =-r s-rrj ; 


sinij) 


2  2  2  W 

k,         t(l  -  V     -  o^,    +  20^,^^'    -  u1    ) 

(1   -   uV    (1    -  u'2)*5 


Gkln(u;k') 


j-  (3u     -   1)1    f        (u')    +  bku   f        (y') 


•    |ck   +  —  (1   -  u   )  J    f    (u*) 


♦  |^««*V)*  (bk-rv)«n«",) 


k 


(3.53) 
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where  the  functions  f  (u ' )  and  g  (u  ' )  are  defined  by 

fn(u')  =  [  u'n  cos-1fa-bu.',]  du1  (3.54) 


and 


with 


gn(u')  =  |  u'n  (c  +  du'  -  V'2)H   du'  ,  (3.55) 


a  =■  ^/(l-u2)4,  (3.56) 

b  =  u/d-u2)*5,  (3.57) 

c  =  1  -  u2  -  u)2,,  (3.58) 

and      d  =  2oV.ll.  (3.59) 

The  integrals  of  Eqs.  (3.54)  and  (3.55),  although  somewhat  complex, 
can  be  evaluated  analytically.   In  the  following  two  subsections  we  will 
examine  the  results  of  these  integrals  for  n  =  0,  1,  and  2.   These 
results  are  sufficient  for  using  linear  or  quadratic  interpolation  for 
the  angular  flux.   If  a  piecewise  linear  representation  of  the 
scattering  cross  section  is  used  (i.e. ,  a,  =  0) ,  these  results  also 
permit  the  use  of  cubic  interpolation  for  the  flux. 

3.3.1  Evaluation  of  the  Functions  f"(u') 

The  argument  of  the  arccosine  in  the  integrand  of  Eq.  (3.54)  is 
singular  at  u'  =  ±1.   From  examination  of  Figs.  3.1  -  3.3,  it  is  clear 

that,  when  u1  =  1,  <J>  =  0  or  tt  (except  when  cu  is  exactly  equal  to  u)  . 

k£n 
Thus  for  u'  =  1,  G    (u;k')  is  given  by  case  3  only  if  u).  ,  =  u.   If  u  , 

=  u,  it  is  clear  from  Eqs.  (3.56)  and  (3.57)  that  a  =  b.   The 

singularity  in  the  integrand  of  Eq.  (3.54)  can  then  be  removed  by 

factorization,  i.e. 
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a  -  bu' |  a(l-u') 


=  0  .         (3.60) 
U'-l  V-l 


(l-u'2)*5!      (i+u')'2(i-u') 


The  singularity  at  U1  =  -1  can  be  removed  in  the  same  manner. 

With  the  singularities  removed,  the  functions  f  (u')  of  Eq.  (3.54) 
can  be  evaluated  analytically  using  integration  by  parts  and  integral 
tables.   '    The  following  results  are  obtained  for  f  (u'): 

f°(u')  =  u'h^u')  +  I  h2(u')  +  h3(u')  -  h4(u'),  (3.61) 

ll (V')   ■!{ll,2h1(u')  +  [at  "3br  )  h2(y')  -h3(y') 


h4(y')  +  ahjOl')}  ,  (3.62) 


f2(V)   =j  jy'^Oi')  +  f^tpOt2? 


2b2  -  a2  +  1)  +  2]  h2(u') 


+  h3(u')  -  h4(u')  +  |  fay'  -  2b  +  3atpi  h^y')},        (3.63) 


f3(u')  =  ^  |w'\(u')  +  *-§•  [5bt2p  -  3bs]  h 


2 

,<U*) 

2r" 


-  -^r   [3t2p  -  3a2  +  2r2  ♦  1]  h  (p')  -  h-(u') 
2r 

-  h4(y')  -  |  [u*  +  3tp]  h5(u') 
+|  [2u'2  +  5tpu'  -  4sp  +  15t2p2  +  6]  h  (u')^  , 

and 
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tV)  =}{u'\(u')  +  -S"  f2^  ~  b)r5t2p  -  3s]h2(y') 


+  -Ar  [3t2p  -  2b2  +  2r2  -  s]  h  (y ' ) 
2r 


-  32.  [3t2p  -  s]  h.(p')  +  h.(u')  -  h.(y') 
8r 

♦  [?atP24~  4bP)  f5—  +  5tu'  "  4s  ♦  15t2p)  h5(y') 


,3 

(M*)V,     (3.65) 


♦  (f-Jpjh,.  ♦3tp]h5(„'>  +  [a^-b)h5 

where 

h.(u')  =  cos'l  a  ~  &C],  (3.66) 

1         l(l-u,2)*J 


h2(p')5.l.-1[-(b2^2';'b],  (3.67) 

2  I    (l-aW)"4  > 


.,,,_!  a+b   ,  -lf(ab+b2+l)(y'+l)  -  (a+b)2"! 
h.(u')  =  -z sin   5 5-r , 

3       2  |a+b|       I   Ip'+ll  (1-a2  +  b2)^  > 


(3.68) 


h.(u')  3  l^^sin-f  ^-b2-l)(p'-l)  -  (a-b)2]   >     (369) 
4       2  |a-b|       I-    lu'-ll  (l-a2+bV    1 


h,(p")  E  /-(b2*Ou''  ♦  2ab"'  -  ^la  ,  (3.70) 

5  b  +  1 


p  E  (b2+l)_1  ,  (3.71) 


52 

r  S  (b^l)15  ,  (3.72) 

s  =   a2  -  1,  (3.73) 

and 

t  =  ab  .  (3.74) 


3.3.2  Evaluation  of  the  Functions  g  ("') 

The  functions  g  (u')  of  Eq.  (3.55)  can  be  evaluated  using  integral 
tables.    The  following  results  are  obtained  for  g  (u'): 

gV).  (2U'-d)(c^--U-2)^di^c  vul)>         (3-75) 


1   ,.  _   (c-t-dy-U'2)3/2   d(d-2u')(c+dU'-U'2)'1 


d(d^4c)  M"'>.  <3-7« 

lb     b 


2,r,  -_[»•♦  "1  (c^'-^'2)3/2  +  Sd^c  0 


16 


and 


where 


i'2/„  .j,ii_  ii  i2%3/2 


g 


3(w.)  ■-"'  C^'-'-  >    +  ^gV)  +  £gV>.   (3.78) 


10 


h.(P')  -  sin"1  -f2^  .  (3.79) 


It  can  be  shown  that  the  function  h.(U')  defined  by  Eq.  (3.79)  is 
equivalent  to  the  function  h„(u')  defined  by  Eq.  (3.67). 
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3.4  Final  Form  of  the  Scattering  Source  Term 

With  the  results  of  sections  3.1  -  3.3,  the  final  form  of  the 
scattering  source  term  can  be  written  as  (cf.  Eqs.  (3.31)  and  (3.47)) 

NM+1 

where  , 

max  M+4  N 
F     (u  -u)  S  2  t        I        I      F^lnOi).  (3. 80 

s  B  k=l  £=1  n=*0 

It  should  be  noted  that  the  form  of  Eq.  (3.80)  is  exactly  the  same  as 
the  form  of  the  exact  kernel  scattering  source  term  (Eq.  (2.48)).   Thus 
any  discrete  ordinates  code  based  on  the  exact  kernel  method  can  be  used 
with  the  above  results. 

In  the  next  chapter,  several  slab  albedo  problems  are  examined. 
Specifically,  results  for  these  problems  obtained  using  the  approximated 
scattering  kernel  method  are  compared  to  those  obtained  using  the  exact 
kernel  and  Legendre  expansion  methods. 
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4.   USE  OF  THE  APPROXIMATED  SCATTERING  KERNEL  METHOD 
IN  DISCRETE-ORDINATES  TRANSPORT  CALCULATIONS 


In  the  previous  chapter,  a  new  method  of  evaluating  the  scattering 
source  term  in  discrete-ordinates  transport  calculations  was  examined. 
The  purpose  of  deriving  this  method  was  to  find  a  means  of  evaluating 
the  scattering  source  term  which,  for  highly  anisotropic  problems, 
provides  better  accuracy  than  does  the  Legendre  expansion  method  and 
which  does  not  suffer  from  the  problem  of  angular  coverage  associated 
with  the  exact  kernel  method. 

In  this  chapter  we  will  apply  this  new  method,  termed  the 
approximated  scattering  kernel  method,  to  several  one-speed  slab  albedo 
transport  problems  and  compare  the  results  obtained  to  those  obtained 
using  the  exact  kernel  and  Legendre  expansion  methods.   However,  before 
examining  the  results  of  these  calculations,  let  us  consider  the  general 
form  of  the  slab  albedo  problem. 

4.1  The  Slab  Albedo  Problem 

Consider  a  source-free,  infinite,  homogeneous  slab  of  thickness  T 

mean-free-path  lengths  in  a  vacuum.   The  slab  is  illuminated  on  one  face 

by  a  monodirectional  beam  of  particles.   This  problem  has  been  widely 

studied  *   '   and  is  known  as  the  "slab  albedo  problem".  Let  the 

incident  flux  at  x  ■  0  be  in  direction  u   >  0,  where  u  is  one  of  the 

s  s 

discrete  ordinates  {u  }.   The  boundary  conditions  are  then 

0(0, u)  =  6(u-ug)  ,  u  >  0  (4.1) 

and 
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♦(t.m)  =0       ,  u  <  0  (4.2) 

where  6(u-u  )  is  Che  Dirac  delta  distribution. 

s 

For  the  simple  case  of  isotropic  scattering,  the  exact  analytical 

solution  for  the  reflected  and  transmitted  angular  fluxes  can  be 

34 
obtained  by  using  the  X-Y  functions  of  radiative  transfer.    In  terms 

of  these  X  and  Y  functions,  the  reflected  and  transmitted  angular  fluxes 

are  given  by 

CPs 
*(°'1J)   "   ofuJ,  •>    [X<W)X(V  )   "  Y(u)Y(u  )]      ,      v  <  0  (4.3) 

2(u+u   )  s  s 

s 

and 

*(T>U>    "   77?.   ,.   ^    rY(u)X(p   )    -  X(u)Y(l»  )]      .      V   >  0    .         (4.4) 
2(u-ii   )  s  s 


Values  of  X(u)  and  Y(y)  are  tabulated  in  Ref.  35  for  various  values  of  c 
and  T. 

For  both  isotropic  and  anisotropic  scattering,  the  diffusely 

reflected  and  transmitted  angular  fluxes  obey  an  important  reciprocity 

34 
principle.     [The  diffusely  reflected  and  transmitted  fluxes  consist 

only  of  those  particles  which  have  undergone  one  or  more  scattering 

processes.   Thus  the  uncollided  transmitted  flux  in  the  source 

direction,  4(0, u  )exp(-x/u  ),  is  not  included  in  the  diffusely 

transmitted  flux  *(t,u  ).]   The  reciprocity  principle  can  be  stated  as 

follows : 

u*(0;u  ;-u)  =  u  *(0;u;-u  )  (4.51 

s        s        s 

and 

u*(t;u  ;u)  =  u  4>(t;u;u  )  (4.6) 
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where  <&(x;u  ;u)  is  the  angular  flux  at  position  x  in  direction  u  due  to 

an  incident  flux  at  x  =  0  in  direction  u  . 

s 

The  difference  between  the  diffusely  transmitted  flux  and  the  total 
transmitted  flux  depends  on  the  source  angle  of  the  incident  flux  as 
well  as  the  slab  thickness.   In  some  cases,  such  as  a  normally  incident 
flux  on  a  thin  slab,  the  uncollided  transmitted  flux  may  be  much  greater 
than  the  diffusely  transmitted  flux.   Since  we  are  primarily  interested 
in  the  effects  of  scattering  source  calculations  in  this  study,  we  will 
compare  only  the  diffuse  fluxes  in  the  problems  of  Sections  4.2  and  4.3. 

4.1.1  Incident-Flux  Normalization  in  Discrete-Ordinates  Transport 
Calculations 

To  obtain  the  proper  angular  fluxes  when  performing  discrete- 

ordinates  transport  calculations,  it  is  necessary  to  normalize  the 

incident  flux  *(0,u  )  in  accordance  with  the  boundary  condition  given  in 
s 

Eq.  (4.1).  The  presence  of  the  delta  distribution  in  Eq.  (4.1)  implies 
that,  for  u  >  0 

4(0, u)  =  0  ,  u  J   u  (4.7) 

s 

and         . 

dy  *(0,u)  -  1  .  (4.8) 

0 

The  specific  form  of  the  source  normalization  is  dependent  on  the 
method  used  to  evaluate  the  integral  of  Eq.  (4.8).   When  numerical 
quadrature  is  used  (as  in  the  exact  kernel  or  Legendre  expansion 
methods),  Eq.  (4.8)  is  evaluated  as 
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NM+1 

I       w  *(0,u.)  =  1  (4.9) 

i-1        x 


where  NM+1  is  the  number  of  discrete  ordinates.   Because  the  incident 

flux  is  zero  in  all  directions  other  than  u  >  the  effective  source 

s 

strength  or  normalization  required  to  satisfy  Eqs.  (4.7)  and  (4.8)  is 
thus  (for  y  >  0) 


«(0,y±)  = 


0    ,   y  l<  u 

x  1    S  (4.10) 


When  the  approximated  scattering  kernel  method  is  used,  the  source 
normalization  is  somewhat  more  complex.   Substitution  of  Eqs.  (3.28)  and 
(3.29)  into  Eq.  (4.8)  yields 


1  NM+1   M   '   jN+1    y-y 


1=1    [  I    \     '   I   — 14  H([y  -y"  ][y   -y  ])  ♦(O.w  )  dy  .(4.11) 
'        1-1  j-1  k=jN-N+l  ui  Mkl 
U  *■   k^i        > 

Since  the  incident  flux  is  zero  in  all  directions  except  y  ,  Eq.  (4.11) 

can  be  simplified  to 

rl       M   I    jN+1    M-V  _  _ 

1  =  ♦(O.u  )     I      I    — S-  H([y  -y  ] [y   -y  ])  dy  .    (4.12) 
s  )q     j_i  k=jN-N+l  ys  Uk      S  J   2    L     S 

Because  the  integrand  of  Eq.  (4.12)  is  zero  outside  the  interval 
[y .  ,y .  .]  for  each  j,  we  can  rewrite  Eq.  (4.12)  as 
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>1 


j-1        J   J        '      k=jN- 


jN+1 


Ma 


du 


(4.13) 


so  that  the  source  normalization  becomes 


KO.iiJ 


I  H([u,-5  ][y  rus 

j=i 


J+11N+ 


jN+1 


n 

N-S 

Ms 


k=jN-N+l  Us~Uk 


du 


-1 


(4.14) 


4.1.2  Treatment  of  the  Discontinuity  in  the  Angular  Flux  at  u=0 
It  may  be  recalled  from  Sec.  2.3.1  that  the  angular  flux  is 
discontinuous  at  u  =  0  at  any  interfaces  in  a  slab.   For  a  homogeneous 
slab,  the  only  interfaces  are  the  slab  surfaces.   The  angular  flux  is 
therefore  continuous  throughout  the  interior  of  the  slab  and 
discontinuous  at  x  =  0  and  x  =  T.   Because  the  angular  flux  is  treated 
as  a  piecewise  continuous  function  in  the  approximated  scattering  kernel 
method,  this  discontinuity  can  be  expected  to  produce  some  error  in  the 
calculation  of  the  reflected  and  transmitted  angular  fluxes. 

One  means  by  which  this  error  can  be  minimized  is  to  place  the 
ordinates  immediately  on  each  side  of  u  ■  0  as  close  to  0  as  is 
practical.   However,  it  must  be  remembered  that  the  number  of  spatial 
nodes  required  to  ensure  positivity  of  the  calculated  angular  fluxes 
increases  as  the  absolute  value  of  the  ordinate  closest  to  u  =  0 
decreases  [cf.  Eq.  2.38)].   This  restriction  places  a  practical 
limitation  on  how  close  to  zero  the  discrete  ordinates  can  be  located. 

An  approach  to  eliminate  the  effects  of  the  discontinuity  at  u  =  0 
is  to  let  the  point  u  =  0  be  one  of  the  discrete  ordinates.   More 
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specifically,  by  letting  u  =  0  be  one  of  the  boundaries  u.  in  the 
piecewise  polynomial  expansion  of  the  angular  flux,  the  scattering 
source  term  will  not  be  integrated  across  u  =  0.   This  approach  is 
similar  to  that  utilized  in  the  DPW  method,  in  that  the  scattering 
source  term  is  evaluated  in  the  two  half-ranges  [-1,0]  and  [0,1]  rather 
than  the  full  range  [-1,1], 

While  this  approach  would  eliminate  integration  across  the  flux 
discontinuity,  it  introduces  another  problem,  namely,  the  evaluation  of 
the  angular  flux  at  y  =  0.   The  standard  discrete-ordinates  equations 
[Eqs.  (2.36)  and  (2.37)]  cannot  be  used  for  u  =  0.   Hence  an  alternative 
method  of  evaluating  the  angular  fluxes  would  have  to  be  utilized  to 
calculate  $(x,0) . 

If  the  flux  around  u  =  0  were  smooth,  one  could  interpolate  the 
flux  at  0  with  a  reasonable  degree  of  accuracy.   When  the  flux  is 
discontinuous,  however,  interpolation  could  produce  very  poor  results. 
Essentially,  then,  interpolation  of  the  flux  at  u  =  0  offers  no 
advantage  over  integrating  across  u  =  0.   As  an  alternative  approach  one 
could  evaluate  the  flux  at  p  =  0+e  and  u  =  0-e,  where  e  is  a  very  small 
number,  by  extrapolating  the  flux  from  the  positive  and  negative 
directions.   Under  certain  conditions,  this  approach  could  work  very 
well.   Consider,  for  example,  the  transmitted  flux  emerging  from  a  slab 
at  x  =  x.   The  angular  flux  in  this  case  is  zero  for  \i   <   0,  and  non-zero 
for  u  >  0.   Extrapolation  of  the  flux  at  0-e  would  produce  the  correct 
value  (i.e.,  *(x,0-e)  =  0).   The  accuracy  of  the  extrapolation  of 
$(t,0+e)»  however,  would  depend  on  the  behavior  of  the  flux  near  U   =  0, 
as  well  as  on  the  location  of  the  discrete  ordinates  used  to  extrapolate 
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the  flux.   If  the  flux  were  smooth  near  u  =  0,  the  extrapolation  could 
be  very  accurate.   However,  if  the  flux  varied  rapidly  near  u  =  0, 
extrapolation  could  yield  a  very  poor  estimate  of  *(x,0+e). 

In  the  transport  calculations  performed  in  this  work,  it  was  found 
that  the  best  results  were  usually  obtained  by  using  the  full-range 
expansion  of  the  scattering  source  term  (and  thus  ignoring  the 
discontinuity  at  the  slab  surfaces).   In  some  problems,  the  use  of  flux 
extrapolation  produced  results  which  were  slightly  more  accurate  than 
those  obtained  without  extrapolation.   However,  in  most  cases  the  use  of 
flux  extrapolation  produced  no  improvement,  and,  in  fact,  often  led  to 
significant  errors.   It  was  therefore  decided  that  the  use  of  flux 
extrapolation  should  be  avoided  because  it  is  not  a  reliable  method  of 
improving  the  accuracy  of  transport  calculations. 

4.2  Results  for  Isotropic  Scattering 

The  first  problem  considered  is  the  slab  albedo  problem  with 

isotropic  scattering.   The  transport  medium  is  a  one  mean-free-path  slab 

with  c  =  1.   Figure  4.1  shows  the  diffusely  reflected  and  transmitted 

angular  fluxes  resulting  from  unit  incident  fluxes  in  directions 

u  =  1.0,  u  =  0.566331,  and  u  =  0.66838,  as  calculated  using  the  X-Y 
s        s  s 

functions. 

The  reflected  and  transmitted  fluxes  for  the  normally  incident  case 
were  first  calculated  using  each  of  the  scattering  source  treatments 
(exact  kernel,  Legendre  expansion,  and  approximated  scattering  kernel) 
with  a  Lobatto-12  discrete-ordinates  set.   The  approximated  scattering 
kernel  calculations  used  an  additional  ordinate  which  was  needed  to 
obtain  an  integral  number  of  flux  intervals  for  piecewise  quadratic  and 
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cubic  flux  expansions.   This  additional  ordinate  was  placed  at  u  =  0.04 
in  an  attempt  to  minimize  the  error  produced  by  the  discontinuity  at 
u  =  0.   The  results  of  all  the  calculations  were  within  about  1%  of  the 
fluxes  calculated  using  the  X-Y  functions.   The  fluxes  calculated  using 
the  approximated  scattering  kernel  method  with  linear,  quadratic,  and 
cubic  flux  expansions  were  actually  somewhat  more  accurate  than  the 
fluxes  calculated  using  the  exact  kernel  and  Legendre  expansion  methods. 
This  difference  may  be  due  at  least  in  part  to  the  extra  ordinate  used 
in  the  approximated  scattering  kernel  calculations. 

An  obvious  means  of  reducing  the  error  in  the  calculated  fluxes  is 
to  increase  the  number  of  discrete  ordinates.   To  illustrate  this 
effect,  the  transport  calculations  were  next  carried  out  for  the 
normally  incident  case  with  a  Lobatto-24  discrete-ordinates  set.   The 
approximated  scattering  kernel  calculations  again  utilized  an  additional 
ordinate  at  u  =  0.04.   For  each  of  the  scattering  source  treatments 
used,  nearly  all  the  angular  fluxes  were  within  0.25%  of  the  X-Y 
results. 

Finally,  calculations  were  performed  with  the  Lobatto-24  set  for 
source  angles  of  0.566331  and  0.066838.   For  the  source  angle  of 
0.566331,  the  angular  fluxes  calculated  using  each  of  the  scattering 
source  treatments  were  again  nearly  all  within  0.25%  of  the  X-Y  fluxes. 
For  the  source  angle  of  0.066838,  the  error  was  slightly  greater,  with 
all  the  fluxes  within  about  0.5%  of  the  X-Y  results. 
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4.3  Results  for  Anisotropic  Scattering 

We  next  consider  three  problems  which  exhibit  increasing  degrees  of 
scattering  anisotropy.   The  transport  medium  in  each  problem  is  again  a 
one  mean-free-path  slab  with  c  =  1. 

Problem  4.3.1 

The  first  anisotropic  scattering  problem  considered  has  a 
scattering  cross  section  whose  range  of  angular  support  is  u  e  [0,  1]. 

This  cross  section,  which  is  shown  in  Fig.  4.2,  is  representative  of  the 

36 
scattering  of  neutrons  by  hydrogen. 

The  reflected  and  transmitted  angular  fluxes  resulting  from  a  unit 
incident  flux  were  calculated  for  source  angles  of  1.0,  0.566331,  and 
0.066838  using  each  of  the  methods  discussed  for  the  evaluation  of  the 
scattering  source  term.   A  Lobatto-24  discrete-ordinates  set  was  used  in 
all  the  calculations.   The  approximated  scattering  kernel  calculations 
with  the  quadratic  and  cubic  flux  expansions  contained  an  additional 
ordinate  which  was  needed  to  obtain  an  integral  number  of  flux 
intervals.   This  additional  ordinate  was  placed  at  u  =  0.04  in  an 
attempt  to  minimize  the  error  incurred  by  interpolating  across  u  =  0. 

The  results  of  the  exact  kernel  calculations,  which  were  judged  to 
be  the  most  accurate  on  the  basis  of  reciprocity,  are  shown  in  Fig.  4.3. 
The  percent  deviations  in  the  angular  fluxes  as  calculated  by  the 
Legendre  expansion  method  and  the  approximated  scattering  kernel  method 
with  linear,  quadratic,  and  cubic  flux  expansions  are  plotted  in  Figs. 

4.4  -  4.6. 

It  can  be  seen  from  Figs.  4.4  -  4.6  that  the  fluxes  calculated  by 
the  approximated  scattering  kernel  method  with  a  linear  flux  expansion 
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give  consistently  good  results.   Nearly  all  of  the  fluxes  calculated  by 
this  method  are  within  1%  of  the  exact  kernel  fluxes. 

Use  of  quadratic  and  cubic  flux  expansion  flux  expansions  in  the 
approximated  scattering  kernel  method  led  to  rather  erratic  results. 
These  methods  often  led  to  fluxes  which  were  badly  over-  or 
underestimated.   The  reason  for  this  behavior  is  discussed  in  Section 
4.5. 

The  fluxes  calculated  using  an  eighth-order  Legendre  expansion  of 
the  scattering  cross  section  exhibited  oscillatory  behavior* 
particularly  for  the  reflected  fluxes.   These  oscillations  were  worst 
for  the  normally  incident  case. 

Problem  4.3.2 

The  second  anisotropic  scattering  problem  considered  has  a 
scattering  cross  section  whose  range  of  angular  support  is  u  e  [0.8,  1]. 
This  fictitious  piecewise  linear  cross  section,  which  is  representative 
of  an  in-group  scattering  cross  section  for  a  multigroup  energy 
structure,  is  shown  in  Fig.  4.7. 

The  reflected  and  transmitted  angular  fluxes  resulting  from  a  unit 
incident  flux  were  calculated  for  source  angles  of  1.0,  0.566331,  and 
0.066838  using  the  same  discrete-ordinates  set  as  in  Problem  4.3.1.   The 
results  of  these  calculations  are  shown  in  Figs.  4.8  -  4.10.   Many  of 
the  fluxes  calculated  using  the  Legendre  expansion  and  approximated 
scattering  kernel  methods  show  significant  deviations  from  the  exact 
kernel  fluxes. 

The  fluxes  calculated  using  the  Legendre  expansion  method  generally 
exhibit  oscillatory  behavior.   The  oscillations  are  most  severe  for  the 
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normally  incident  flux,  in  which  case  most  of  the  reflected  fluxes  were 
actually  calculated  to  be  negative.   For  the  other  two  source  angles, 
the  oscillations  were  again  most  pronounced  for  the  reflected  fluxes. 
Note  that  these  results  are  qualitatively  similar  to  those  obtained  in 
Problem  4.3.1  (cf.  Figs.  4.4  -  4.6). 

The  fluxes  calculated  using  the  approximated  scattering  kernel 
method  also  exhibit  the  same  type  of  behavior  shown  in  the  results  of 
Problem  4.3.1.   Calculations  using  a  linear  expansion  of  the  flux 
generally  yielded  the  closest  agreement  with  the  exact  kernel  fluxes. 
The  deviation  between  the  two,  however,  was  much  more  significant  than 
in  Problem  4.3.1.   For  the  normally  incident  case,  the  reflected  fluxes 
calculated  using  a  linear  flux  expansion  were  up  to  40%  higher  than  the 
exact  kernel  fluxes.   The  transmitted  fluxes  were  up  to  20%  higher  than 
those  calculated  using  the  exact  kernel  method.   For  the  other  two 
source  angles,  the  reflected  fluxes  calculated  using  a  linear  flux 
expansion  were  generally  within  107.   of  the  exact  kernel  fluxes,  and  the 
transmitted  fluxes  within  5%. 

The  fluxes  calculated  using  quadratic  and  cubic  flux  expansions 
again  yielded  rather  erratic  results.   It  can  be  seen  from  Fig.  4.8  that 
the  quadratic  and  cubic  flux  expansions  produced  unrealistic  transmitted 
flux  distributions  in  the  normally  incident  case.   The  oscillatory 
behavior  of  these  fluxes  is  most  likely  due  to  the  inability  of 
piecewise  quadratic  and  cubic  expansions  to  model  the  angular  flux 
distributions  realistically.   This  phenomenon  is  discussed  more 
thoroughly  in  Sec.  4.5. 
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For  the  incident  flux  at  0.566331,  the  quadratic  expansion  results 
were  in  better  agreement  with  the  exact  kernel  fluxes  than  were  those 
obtained  using  any  of  the  other  methods  tested.   The  cubic  flux 
expansion  yielded  the  poorest  results  for  this  case.   Finally,  the 
fluxes  calculated  using  quadratic  and  cubic  flux  expansions  gave  rather 
poor  results  for  the  source  angle  of  0.066838.   Calculations  using  the 
cubic  flux  expansion  actually  predicted  negative  angular  fluxes  near 
II  -  1. 

Problem  4.3.3 

The  final  problem  considered  has  scattering  that  is  even  more 
anisotropic  than  the  previous  problem,  with  a  scattering  cross  section 
whose  range  of  angular  support  is  to  e  [0.95,  1],   This  fictitious 
piecewise  linear  scattering  cross  section,  which  is  shown  in  Fig.  4.11, 
is  typical  of  in-group  neutron  scattering  from  light  elements  or  with  a 
fine-energy-group  structure  (cf.  Fig.  2.3). 

The  reflected  and  transmitted  angular  fluxes  resulting  from  a  unit 
incident  flux  were  calculated  for  source  angles  of  1.0,  0.544182,  and 
0.044247.  A  Lobatto-36  discrete-ordinates  set  was  used  in  all  the 
calculations.   The  approximated  scattering  kernel  calculations  with 
quadratic  and  cubic  flux  expansions  again  contained  an  additional 
ordinate  at  u  =  0.04. 

Figure  4.12  shows  the  reflected  and  transmitted  angular  fluxes  for 
the  normally  incident  case  as  calculated  using  the  exact  kernel  method 
and  the  approximated  scattering  kernel  method  with  a  linear  flux 
expansion.   Calculations  using  the  Legendre  expansion  method  and  the 
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approximated  scattering  kernel  method  with  quadratic  and  cubic  flux 
expansions  yielded  oscillatory  and  negative  flux  values  for  both  the 
reflected  and  transmitted  fluxes. 

Note  that  the  fluxes  calculated  using  the  linear  flux  expansion 
significantly  overpredict  the  exact  kernel  fluxes  over  the  entire 
reflected  range  and  much  of  the  transmitted  range.   This  behavior  is 

qualitatively  similar  to  that  shown  in  the  results  of  Problem  4.3.2  (cf. 

20 
Fig.  4.8).   It  is  of  interest  to  note  that  Odom   found  that  small 

inaccuracies  in  scattering  source  calculations  can  produce  relatively 

large  inaccuracies  in  the  angular  flux,  particularly  in  reflected 

fluxes. 

The  reflected  and  transmitted  angular  fluxes  as  calculated  by  the 
exact  kernel  method  and  the  approximated  scattering  kernel  method  with 
linear,  quadratic,  and  cubic  flux  expansions  are  shown  in  Fig.  4.13  for 
the  incident  flux  in  direction  0.544182.   The  results  obtained  using 
linear  and  quadratic  expansions  of  the  flux  are  in  good  agreement  with 
the  exact  kernel  fluences,  particularly  for  the  transmitted  fluxes.   The 
fluxes  calculated  using  a  cubic  expansion  show  poor  behavior  in  the 
transmitted  directions.   Calculations  using  the  Legendre  expansion 
method  again  produced  oscillatory  and  negative  angular  fluxes. 

The  reflected  and  transmitted  angular  fluxes  as  calculated  by  the 
exact  kernel  method  and  the  approximated  scattering  kernel  method  with 
linear  and  quadratic  flux  expansions  are  shown  in  Fig.  4.14  for  the 
incident  flux  in  direction  0.044247.   Both  of  the  approximated 
scattering  kernel  flux  profiles  are  in  good  agreement  with  the  exact 
kernel  fluxes.   Calculations  using  the  Legendre  expansion  method  and  the 
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approximated  scattering  kernel  method  with  a  cubic  flux  expansion 
yielded  oscillatory  and  negative  flux  values. 

4.4  Angular  Redistribution  in  the  Approximated  Scattering  Kernel  Method 

It  can  be  seen  from  the  results  of  the  problems  in  Sec.  4.3  that 
the  approximated  scattering  kernel  method  provides  better  accuracy  than 
the  conventional  Legendre  expansion  method  for  highly  anisotropic 
problems.   This  increased  accuracy  was  the  first  objective  in  developing 
the  approximated  scattering  kernel  method.   In  this  section  we  will 
examine  the  second  objective,  namely,  the  elimination  of  the  angular 
coverage  problem  which  limits  the  exact  kernel  method  for  low-order 
quadrature. 

The  failure  of  the  exact  kernel  method  to  produce  angular 
redistribution  of  scattered  particles  for  coarse  quadrature  sets  is  due 
to  the  discretization  of  the  angular  variable  u.   The  approximated 
scattering  kernel  method,  on  the  other  hand,  treats  the  incident 
particle  direction  u'  as  a  continuous  variable  [cf.  Eqs.  (3.47)  and 
(3.48)].   Thus,  instead  of  only  those  particles  in  some  discrete 
direction  u.  contributing  to  the  scattering  source  term  in  some  other 
direction  u .,  all  those  particles  with  incident  direction  u'  in  some 
specified  range  are  considered.   This  specified  range  is  determined  by 
the  breakpoints  described  in  Sec.  3.2.2. 

The  ability  of  the  approximated  scattering  kernel  method  to  allow 
angular  redistribution  for  problems  in  which  the  exact  kernel  method 
fails  to  can  be  illustrated  by  the  following  example.   Consider  a  unit 
flux  normally  incident  on  the  slab  described  in  Problem  4.3.3.   Let  the 
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discrete  ordinates  be  given  by  a  LobaCto-10  quadrature  set.   It  can  be 
shown  by  use  of  Eqs.  (2.51)  -  (2.53)  that  the  Lobatto-10  quadrature  set 
has  "zero  order  of  angular  coverage"  and  hence  is  unsuitable  for  exact 
kernel  transport  calculations.   The  normally  incident  particles  fail  to 
scatter  into  any  other  discrete  directions  and  hence  the  fluxes  in  all 
other  directions  are  calculated  to  be  zero  (see  Table  4.1).   Use  of  the 
approximated  scattering  kernel  method,  on  the  other  hand,  allows  angular 
redistribution  from  the  source  direction  to  all  other  directions  (see 
Table  4.1). 

Although  use  of  the  approximated  scattering  kernel  method  will 
allow  angular  redistribution  with  any  discrete-ordinates  set,  the 
accuracy  obtained  may  be  poor.   For  example,  the  fluxes  calculated  using 
the  approximated  scattering  kernel  method  in  the  example  above  show 
significant  difference  from  those  calculated  using  the  exact  kernel 
method  with  a  Lobatto-36  quadrature  set.   The  reasons  for  this 
discrepancy  are  discussed  in  the  next  section. 

4.5  Analysis  of  the  Approximated  Scattering  Kernel  Results 

It  can  be  seen  from  the  problems  examined  in  Sections  4.2  -  4.4 
that  the  use  of  the  approximated  scattering  kernel  method  in  discrete- 
ordinates  transport  calculations  sometimes  leads  to  very  accurate 
results,  and  yet  at  other  times  results  in  very  poor  estimates  of 
angular  fluxes.   In  this  section,  we  will  examine  the  factors  which 
affect  the  accuracy  of  the  approximated  scattering  kernel  method. 

The  first,  and  probably  most  important,  factor  to  consider  is  the 
degree  of  anisotropy  in  the  angular  flux.   Examination  of  the  results  of 
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Table  4.1.   The  reflected  and  transmitted  angular  fluxes 
resulting  from  a  unit  flux  normally  incident 
on  a  one  mean-free-path  slab  with  c  =  1  and 
the  scattering  cross  section  of  Fig.  4.11. 
The  discrete  ordinates  are  given  by  a  Lobatto-10 
quadrature  set. 


Cosine  of 

Approximated 

Polar  Angle 

Exact  Kernel 

Scattering  Kernel 

Reference 

(u) 

Flux 

Flux 

Flux3 

-1.0 

0.000 

5.923(-8)b 

1 .333C-13) 

-0.919534 

0.000 

1.184(-6) 

1.602(-11) 

-0.738774 

0.000 

2.107(-5) 

1.657(-9) 

-0.477925 

0.000 

1.916(-4) 

6.166(-8) 

-0.165279 

0.000 

4.813(-4) 

2.433(-7) 

0.165279 

0.000 

7.899(-3) 

8.550(-6) 

0.477925 

0.000 

4.062(-2) 

2.315(-4) 

0.738774 

0.000 

3.809(-l) 

1.518(-2) 

0.919534 

0.000 

3.570 

9.714(-1) 

1.0 

2.371(1) 

1.144(1) 

2.108(1) 

The  reference  flux  values  were  calculated  by  performing  cubic 
spline  interpolation  on  the  exact  kernel  fluxes  calculated  with 
a  Lobatto-36  quadrature  set  (cf.  Fig.  4.12). 


Read  as  5.923  x  10 


84 

the  problems  in  Sections  4.2  and  4.3  reveals  that,  as  the  degree  of 
anisotropy  in  the  angular  flux  increases,  the  difference  between  the 
approximated  scattering  kernel  results  and  the  exact  kernel  results 
becomes  greater.   This  increased  error  is  due  to  the  inability  of 
piecewise  polynomial  interpolation  to  accurately  model  highly 
anisotropic  angular  fluxes.   This  can  be  shown  more  precisely  by 

considering  the  error  associated  with  polynomial  interpolation,  which  is 

,  37 
given  by 


,(n+l).  .  j+n 


where  E(u)  is  the  error  in  interpolating  the  function  at  u ,  f(u)  is  the 
continuous  function  which  is  approximated  by  the  interpolating 
polynomial,  the  u  are  the  n+1  values  of  u  at  which  the  function  is 
known,  and  5  is  a  function  of  u  and  is  within  the  range  of  the  u.,  but 
is  otherwise  unknown. 

Equation  (4.15)  does  not  allow  us  to  calculate  the  interpolation 
error,  for  we  do  not  know  the  function  f(u)  which  describes  the  angular 
flux,  nor  do  we  know  the  value  of  £.   However,  Eq.  (4.15)  does  indicate 
that  the  interpolation  error  is  proportional  to  the  (n+l)'th  derivative 
of  f(u).   This  derivative  would  be  expected  to  increase  as  the  angular 
flux  becomes  more  anisotropic. 

Another  factor  which  affects  the  accuracy  obtained  in  discrete- 
ordinates  calculations  which  utilizes  the  approximated  scattering  kernel 
method  is  the  order  of  flux  expansion  used.   Use  of  a  piecewise  linear 
flux  expansion  yields  results  which  are  often  very  accurate,  and  are 
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always  at  least  physically  realistic.   Use  of  quadratic  and  cubic  flux 
expansions,  on  the  other  hand,  can  lead  to  erratic  results.   In  some 
cases  these  expansions  perform  very  well,  while  at  other  times  they  lead 
to  estimates  of  oscillatory  and  even  negative  angular  fluxes.   This 
behavior  should  not  be  surprising,  as  one  of  the  difficulties  associated 

with  polynomial  interpolation  is  the  oscillatory  behavior  which  it  is 

37 
possible  for  high  order  interpolating  polynomials  to  assume. 

For  example,  consider  the  following  transmitted  fluxes  for  the 

normally  incident  case  in  Problem  4.3.3  as  calculated  by  the  exact 

kernel  method: 


v± 

t(u±) 

0.851155 

0.17756 

0.894266 

0.56950 

0.930378 

1.4096 

0.959209 

5.4667 

0.980532 

13.299 

0.994179 

18.654 

1.0 

252.84 

These  points  are  plotted  in  Fig.  4.15,  along  with  the  piecewise  linear, 
quadratic,  and  cubic  polynomials  that  pass  through  the  exact  kernel 
results.   (The  flux  in  direction  y.  ■  1.0  consists  of  a  diffusely 
transmitted  component  of  magnitude  21.08  and  an  uncollided  component  of 
magnitude  231.76.   Only  the  diffusely  transmitted  flux  is  plotted  in 
Fig.  4.12.) 

Between  a±   =  0.851155  and  u  =  0.959209,  the  flux  exhibits  very 
smooth  behavior  and  all  the  piecewise  fits  produce  essentially  the  same 
results.   Between  v±   =  0.959209  and  u.  =  1.0,  the  results  are  much 
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different.   The  "jump"  in  the  angular  flux  at  u  =  1.0  produces  serious 
oscillations  in  the  quadratic  and  cubic  expansions  for  the  segments 
which  contain  u.  =  1.0.   The  piecewise  linear  expansion  appears  to  be  a 
much  more  realistic  fit  in  this  region. 

Although  there  is  no  way  to  determine  a  priori   when  quadratic  and 
cubic  flux  expansions  will  fail  to  produce  accurate  flux  profiles,  it  is 
obvious  from  the  results  of  Problems  4.3.1  -  4.3.3  (as  well  as  the  above 
example)  that  quadratic  and  cubic  expansions  are  very  sensitive  to  the 
degree  of  anisotropy  in  the  flux.   Thus  one  should  be  very  cautious 
about  using  quadratic  or  cubic  flux  expansions  when  dealing  with 
problems  that  are  likely  to  yield  highly  anisotropic  angular  fluxes. 
The  results  of  Problems  4.3.2  and  4.3.3  show  that  cubic  flux  expansions 
are  particularly  prone  to  error. 

Another  factor  which  affects  the  accuracy  which  can  be  obtained 
when  the  approximated  scattering  kernel  method  is  used  is  the  number  of 
discrete  ordinates  which  are  used.   In  Sec.  4.2,  increasing  the  number 
of  discrete  ordinates  resulted  in  a  decrease  in  error  no  matter  which 
method  was  used  to  evaluate  the  scattering  source.   Similar  comparisons 
(which  are  not  presented  here)  for  the  anisotropic  scattering  problems 
of  Section  4.3  revealed  the  same  behavior.   On  the  other  hand, 
decreasing  the  number  of  discrete  ordinates  in  Problem  4.3.3  led  to  a 

significant  increase  in  error  (cf.  Sec.  4.4).   It  is  of  interest  to  note 

38 

that  Atkinson   states  that,  when  using  piecewise  polynomial 

interpolation,  it  is  often  more  advantageous  to  increase  the  number  of 
interpolation  nodes  then  to  increase  the  degree  of  the  interpolating 
polynomial,  due  to  the  tendency  of  higher  order  polynomials  to  cause 
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large  oscillations  between  the  nodes.   Such  appears  to  be  the  best 
procedure  for  increasing  the  accuracy  of  the  approximated  scattering 
kernel  results. 

An  obvious  means  of  avoiding  the  spurious  oscillations  in  the 
interpolated  flux  profiles  would  be  to  utilize  a  "smoother"  interplation 
technique,  such  as  cubic  spline  interpolation.   However,  the  use  of  a 
technique  such  as  cubic  spline  interpolation  would  require  a  significant 
modification  of  the  approximated  scattering  kernel  method. 

When  Lagrange  interpolation  is  used,  the  angular  flux  in  the 

segment  [5  »U.+,]  is  expanded  as  Tcf.  Eqs.  (3.26)  and  (3.45)] 

jN+1  i  N 

III1  5j  <  V  <  v.+1   . 


1N  +  1  f    N  .  I 

i    {i  cMvv 

t-jH-N+l      4i=0  '     s 


Cubic  spline  interpolation,  on  the  other  hand,  leads  to  an  expansion  of 

the  form 

3 
*  00  -  I  A"!  u"   ,   u.  <  u  <  u    .  (4.17) 

8      n=0  J        3   l 

Comparison  of  Eqs.  (4.16)  and  (4.17)  reveals  why  cubic  spline 

interpolation  cannot  be  used  in  the  approximated  scattering  kernel 

method  as  derived  in  Chapter  3.   In  Eq.  (4.16),  the  quantities  which 

change  with  each  "inner  iteration"  of  the  discrete-ordinates  solution 

algorithm  are  the  angular  flux  values  *  (tO.   The  expansion 

g  i 

coefficients  C  are  dependent  only  on  the  discrete-ordinates  set  and 
hence  are  independent  of  the  iterative  process.   Thus  the  approximated 
scattering  kernel  transfer  matrix  needs  to  be  calculated  only  once  for  a 
particular  transport  medium  and  discrete-ordinates  set. 
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On  Che  other  hand,  if  cubic  spline  interpolation  were  used  [i.e., 
Eq.  (4.17)],  the  quantities  which  change  with  each  flux  iteration  would 
be  the  expansion  coefficients  A-' .   These  coefficients  are  dependent  on 
the  discrete  ordinates  used  as  well  as  on  the  angular  flux  values. 
Hence  the  approximated  scattering  kernel  transfer  matrix  would  have  to 
be  recalculated  with  each  inner  iteration.   While  such  an  approach  is 
possible,  it  would  add  considerably  to  the  computational  effort  of  the 
transport  calculations,  and  has  not  been  pursued  in  this  study. 
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5.   CONCLUSIONS 

In  this  work  a  new  method  of  evaluating  the  scattering  source  term 
in  discrete-ordinates  transport  calculations  has  been  studied.   This 
method,  termed  the  approximated  scattering  kernel  method,  was  developed 
especially  for  use  in  highly  anisotropic  problems  (i.e.,  those  problems 
which  are  characterized  by  highly  anisotropic  angular  fluxes  as  well  as 
by  highly  anisotropic  scattering).  The  objectives  in  developing  the 
approximated  scattering  kernel  method  were  twofold  —  to  achieve  better 
accuracy  than  that  afforded  by  the  conventional  Legendre  expansion 
method,  and  to  eliminate  the  problem  of  angular  redistribution  that 
limits  the  exact  kernel  method  for  low  quadrature  orders. 

The  first  objective,  achieving  better  accuracy  than  the  Legendre 
expansion  method,  was  easily  met  when  the  piecewise  linear  flux 
expansion  model  was  utilized.   For  the  problems  considered  in  Chapter  4, 
use  of  the  approximated  scattering  kernel  method  with  a  piecewise  linear 
flux  expansion  consistently  yielded  angular  fluxes  which  were  as 
accurate  as,  or  more  accurate  than,  the  Legendre  expansion  results.   For 
highly  anisotropic  slab  albedo  problems,  the  approximated  scattering 
kernel  results  were  far  superior  to  the  Legendre  expansion  results, 
which  exhibited  oscillatory  and  even  negative  angular  fluxes. 

The  second  objective  of  the  approximated  scattering  kernel  method, 
the  elimination  of  angular  coverage  problems,  was  also  met.   In  the 
final  example  problem  in  Chapter  4,  use  of  the  approximated  scattering 
kernel  method  with  a  piecewise  linear  flux  expansion  produced  nonzero 
fluxes  while  the  exact  kernel  method  failed  to  achieve  any  angular 
redistribution  at  all.   The  ability  of  the  approximated  scattering 
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kernel  method  to  achieve  angular  redistribution  with  any  number  of 
discrete  ordinates  is  due  to  the  piecewise  polynomial  expansion  of  the 
angular  flux.   In  this  sense,  the  approximated  scattering  kernel  method 
can  be  considered  a  combination  modal  and  nodal  method,  In  that  it  uses 
a  modal  approach  (expansion  of  the  angular  flux  and  the  scattering  cross 
section)  to  generated  nodal  quantities  (the  elements  of  the  approximated 
scattering  kernel  transfer  matrix). 

As  is  the  case  with  both  the  Legendre  expansion  method  and  the 
exact  kernel  method,  the  approximated  scattering  kernel  method  is  not 
without  its  limitations.   Like  the  exact  kernel  method,  it  is  somewhat 
more  cumbersome  to  use  than  the  Legendre  expansion  method,  due  to  the 
necessity  of  generating  and  storing  a  complete  transfer  matrix  rather 
than  just  a  few  expansion  coefficients.   Furthermore,  it  yields  results 
which  are  generally  not  as  accurate  as  the  exact  kernel  method  for  the 
same  discrete-ordinates  set.   In  addition,  the  use  of  piecewise 
quadratic  and  piecewise  cubic  flux  expansions  in  the  approximated 
scattering  kernel  method  leads  to  inconsistent  results.   Although  these 
methods  sometimes  work  well,  they  produce  oscillatory  and  even  negative 
angular  fluxes  in  other  instances.   The  reason  for  these  erratic  results 
is  the  oscillatory  behavior  which  the  piecewise  quadratic  and  cubic 
interpolating  polynomials  can  assume.   It  was  found  that  this 
oscillatory  behavior  is  very  sensitive  to  the  degree  of  anisotropy  in 
the  angular  flux. 

In  general,  it  can  be  concluded  that  the  adequacy  of  the 
approximated  scattering  kernel  method  is  entirely  dependent  on  the 
degree  to  which  the  angular  flux  exhibits  low-order,  piecewise 
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polynomial  behavior.   If  the  angular  flux  can  be  well  modeled  by  a 
piecewise  polynomial  expansion,  the  approximated  scattering  kernel 
method  can  be  expected  to  work  well.   If  the  angular  flux  exhibits 
"non-polynomial"  behavior  (e.g.,  a  highly  anisotropic  profile),  the 
approximated  scattering  kernel  method  may  produce  very  poor  results. 
Thus,  any  factors  which  affect  the  anisotropy  of  the  angular  flux  can 
significantly  affect  the  accuracy  of  discrete-ordinates  transport 
calculations  which  utilize  the  approximated  scattering  kernel  method. 
For  example,  use  of  the  approximated  scattering  kernel  method  with  a 
linear  flux  expansion  in  Problem  4.3.3  yielded  accurate  fluxes  for  a 
grazing  angle  of  incidence  (y  =  0.044247),  yet  greatly  overestimated 
many  of  the  angular  fluxes  for  the  same  problem  when  the  slab  was 
normally  illuminated  (u  =1.0). 

In  spite  of  its  shortcomings,  however,  the  approximated  scattering 
kernel  method  is  a  practical  method  for  scattering  source  term 
calculations  in  discrete-ordinates  transport  calculations.   Provided  a 
piecewise  linear  flux  expansion  is  utilized,  the  approximated  scattering 
kernel  method  yields  positive  (and  apparently  physically  realistic) 
angular  fluxes  for  problems  in  which  the  Legendre  expansion  method 
fails.   Unlike  the  exact  kernel  method,  it  produces  angular 
redistribution  of  scattered  particles  with  any  number  of  discrete 
ordinates.   Furthermore,  by  obviating  the  need  for  numerical  quadrature 
in  the  scattering  source  term  calculation,  it  allows  one  to  choose  a 
discrete-ordinates  set  without  the  usual  constraints  involved  in  that 
process. 
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Further  studies  in  this  area  could  deal  with  ways  of  minimizing  the 
error  associated  with  interpolation  of  the  angular  flux.   One 
possibility  would  be  to  attempt  to  optimize  the  location  of  the  discrete 
ordinates.   This  problem  has  been  widely  studied  for  global  polynomial 
interpolation,  in  which  case  choosing  the  N  discrete  ordinates  as  the 

zeroes  of  the  N'th  order  Chebyshev  polynomial  tends  to  (but  does  not 

37 
always)  minimize  the  maximum  value  of  the  interpolation  error. 

However,  this  procedure  does  not  apply  to  piecewise  polynomial 

interpolation. 

If  one  knew  a  priori   approximately  what  the  angular  flux  profile 
for  a  particular  problem  was,  one  could  minimize  the  interpolation  error 
by  locating  many  discrete  ordinates  where  the  angular  flux  varies 
rapidly.   However,  while  it  may  be  possible  to  choose  such  an  optimal 
discrete-ordinates  set  for  a  particular  problem  (e.g.,  a  flux  incident 
at  some  angle  u  on  a  slab  of  a  certain  thickness  with  a  certain 
scattering  cross  section) ,  a  change  in  any  one  of  the  problem  parameters 
(such  as  the  source  angle)  could  result  in  a  flux  profile  substantially 
different  from  the  original  one,  so  that  the  discrete  ordinates  would  no 
longer  be  optimally  located.   It  should  be  obvious  that  it  would  be 
extremely  inconvenient  as  well  as  inefficient  to  use  a  different 
discrete-ordinates  set  (and  hence  have  to  generate  a  new  approximated 
scattering  kernel  transfer  matrix)  for  every  transport  problem  one 
wished  to  solve. 

Perhaps  a  more  feasible  means  of  improving  the  accuracy  of  the 
approximated  scattering  kernel  method  would  be  to  investigate  other 
functional  representations  for  the  angular  flux.   Such  representations 
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might  employ  exponential  or  logarithmic  terms,  or  even  rational 
functions.   Use  of  such  non-polynomial  representations  would  require 
significant  changes  in  the  approximated  scattering  kernel  method  as 
derived  in  Chapter  3. 

Finally,  the  approximated  scattering  kernel  method  could  be  applied 
to  multigroup  problems  in  plane  or  spherical  geometry.   This  would 
require  no  change  in  the  calculation  of  the  approximated  scattering 
kernel  transfer  matrices,  and  only  minor  changes  in  the  transport  code 
utilized  in  this  work. 
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8.   APPENDIX 
Computer  Codes 

The  three  main  computer  programs  which  were  utilized  in  this  work 
are  listed  in  this  appendix.  A  brief  description  of  each  code  is  given 
here.   The  code  listings  are  supplied  with  comment  cards  to  facilitate 
the  understanding  and  use  of  the  programs.   All  the  programs  are  written 
in  FORTRAN-77.   They  were  run  on  the  Kansas  State  University  Computing 
Center's  NAS  6630,  which  is  functionally  equivalent  to  an  IBM  4831. 

The  first  code  listed  is  EXKERNEL,  which  is  used  to  generate  exact 
kernel  transfer  matrices.   EXKERNEL  uses  a  piecewise  linear  or  piecewise 
quadratic  approximation  of  the  scattering  cross  section.   This  method 
allows  the  exact  kernel  cross  sections  to  be  evaluated  without  recourse 
to  numerical  quadrature.   The  only  input  data  required  are  the 
breakpoints  in  the  scattering  cross  section,  the  values  of  the  cross 
section  at  the  breakpoints  (plus  the  cross  section  values  at  an 
additional  point  in  each  subdomain  if  a  piecewise  quadratic  cross 
section  approximation  is  used) ,  and  a  set  of  polar  quadrature  ordinates 
(i.e.,  the  discrete  ordinates  {u  }) . 

The  second  code  listed  is  ASKERNEL,  which  is  used  to  generate 
approximated  scattering  kernel  transfer  matrices.   For  transfer  between 
any  two  directions  of  a  discrete-ordinates  set,  ASKERNEL  calculates  the 
various  breakpoints  in  the  integration  range  of  the  scattering  source, 
determines  the  form  of  the  integrand  in  each  subrange,  performs  the 
required  integrations,  and  sums  the  various  components  to  obtain  the 
approximated  scattering  kernel.   The  input  requirements  are  basically 
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the  same  as  those  for  EXKERNEL.   In  addition,  the  user  must  specify  what 
order  of  piecewise  flux  expansion  he  wishes  to  use  (N)  as  well  as  the 
number  of  flux  intervals  (M) .   The  sum  NM+1  must  equal  the  number  of 
discrete  ordinates  which  are  used. 

The  final  code  listed  is  TRANS,  which  performs  one-dimensional, 
azimuthally  symmetric  plane  geometry  discrete-ordinates  transport 
calculations.   The  version  of  TRANS  listed  here  uses  the  approximated 
scattering  kernel  technique  for  the  evaluation  of  the  scattering  source 
term.   Modified  versions  of  TRANS  were  used  to  perform  transport 
calculations  with  the  Legendre  expansion  and  exact  kernel  methods. 
The  input  requirements,  which  vary  slightly  according  to  the  version 
used,  are  described  in  the  program  listing. 


EXKERNEL.FOR  10° 

06-25-1986 

C 

C  ***  PROGRAM  NAME:   EXKERNEL  *** 

C 

C     PROGRAM  TO  INTEGRATE  SCATTERING  CROSS  SECTIONS  OVER  THE 

C     AZIMUTHAL  ANGLE  TO  GENERATE  AZIMUTHALLY  INDEPENDENT  EXACT 

C     KERNEL  CROSS  SECTIONS.   A  PIECEWISE  LINEAR  OR  PIECEWISE 

C     QUADRATIC  CROSS  SECTION  EXPANSION  IS  USED  SO  THAT  THE  EXACT 

C     KERNEL  CROSS  SECTIONS  CAN  BE  INTEGRATED  ANALYTICALLY.   THIS 

C     ALLOWS  A  MUCH  MORE  RAPID  COMPUTATION  THAN  DOES  THE  USE  OF 

C     NUMERICAL  QUADRATURE. 

C 

C     MAJOR  VARIABLES  ARE  AS  FOLLOWS: 

C     'NR1         =  NUMBER  OF  REFLECTED  POLAR  QUADRATURE  DIRECTIONS 

C     -NT1         =  NUMBER  OF  TRANSMITTED  POLAR  QUADRATURE  DIRECTIONS 

C     •KMAX'       =  NUMBER  OF  DISTINCT  BREAKPOINTS  IN  THE 

C  DIFFERENTIAL  SCATTERING  CROSS  SECTION  (I.E.,  DO 

C  NOT  COUNT  DUPLICATE  BREAKPOINT  VALUES  MORE  THAN 

C  ONCE) 

C     'W(I)"       =  BREAKPOINTS  IN  THE  ANGULAR  CROSS  SECTION  (PLUS 

C  THE  MIDPOINTS  OF  EACH  CROSS  SECTION  SUBRANGE  IF 

C  A  QUADRATIC  CROSS  SECTION  EXPANSION  IS  USED) 

C     •SIG(I)'     =  CROSS  SECTION  VALUES  AT  THE  W(I)  VALUES 

C     'NORD-       =  ORDER  OF  CROSS  SECTION  FIT.  INPUT  'I'  FOR  A 

C  LINEAR  CROSS  SECTION  FIT  AND  '2'  FOR  A 

C  QUADRATIC  CROSS  SECTION  FIT. 

C     •U(N)'       =  POLAR  QUADRATURE  ORDINATES 

C 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A(64,64) 

C0MM0N/BL0CK1/  AA(3) ,BB(3) ,CC(3) ,W(7) ,SIG(7) 

COMM0N/BL0CK2/  WB(4) ,SIGB(4) ,U(64) 
C 

READ(5.800)  KMAX. NORD 

800  F0RMAT(2(I5)) 
KMAX=KMAX-1 
KMAX1=KMAX*N0RD+1 

C 

C     READ  IN  ANGULAR  BREAKPOINTS  AND  THEIR  CORRESPONDING  SCATTERING 

C     CROSS  SECTIONS 

READ(5,801)  (W(I),I=1.KMAX1) 

READ(5,801)  (SIG(I),I=1,KMAX1) 

801  F0RMAT(4(D18.8)) 
C 

KMAX2=KMAX+1 

DO  1  I=1,KMAX2 

WB(I)=W((I-1)*N0RD+1) 

SIGB(I)=SIG((I-1)*N0RD+1) 
1  CONTINUE 
C 
C     READ  QUADRATURE  ORDINATES 

READ(5,802)  NTOT 

802  F0RMAT(I5) 

READ(5,803)  (U(I),  1=1. NTOT) 

803  F0RMAT(3(D24.15)) 


EXKERNEL.FOR  101 

06-25-1986 

C 

C     CHECK  FOR  SYMMETRY  OF  ANGULAR  MESH.   IF  IT  IS  SYMMETRIC,  ONLY 
C     ONE-HALF  THE  EXACT  KERNEL  CROSS  SECTIONS  NEED  BE  CALCULATED. 
IFLAG=1 
DO  2  I=l,NTOT 

IF(U(I).NE.-U(NTOT+l-I))  THEN 
IFLAG=0 
GO  TO  3 
END  IF 

2  CONTINUE 

3  CONTINUE 
C 

IF(IFLAG.EQ.l)  THEN 

NT0T1=IFIX(  (NTOT+1  )/2.  ) 

ELSE 

NTOT1=NTOT 

END  IF 
C 

C     CALL  SUBROUTINE  TO  EVALUATE  EXPANSION  COEFFICIENTS  FOR  CROSS 
C     SECTION  EVALUATION. 
C 

CALL  COEFF(KMAX,NORD) 
C 

DO  5  I=1,NT0T1 

DO  5  J=l,NTOT 

0MEGA1=U( I)*U( J)+DSQRT( ( 1 .  -U( I )**2)*(  I .  -U(  J}**2) ) 

0MEGA2=U( I)*U(  J)-DSQRT(  { 1 .  -U( I )**2)*(  1 . -U( J)**2) ) 
C 

C     IF  SCATTERING  RANGES  IMPOSED  BY  KINEMATICS  AND  BY  QUADRATURE 
C     SET  DO  NOT  OVERLAP,  CROSS-SECTION  IS  ZERO 
C 

IF(0MEGA1 .  LT .  W(  1 )  .  OR .  0MEGA2 .  GT .  W(KMAX2)  )  THEN 
A(I,J)=0. 
GO  TO  4 

ELSE 

A(I,J)=FNSIGG(I,J,KMAX) 

END  IF 
C 

4  CONTINUE 

IF(IFLAG.EQ.l)  A(NTOT+1-I,NTOT+1-J)=A(I, J) 
C 

5  CONTINUE 

WRITE(6,900)  ((A(I,J),J=1,NT0T),I=1,NT0T) 
900  F0RMAT(4(D18.8)) 

STOP 

END 
C 
C 
C 
C 

C     SUBROUTINE  TO  EVALUATE  EXPANSION  COEFFICIENTS  FOR  LINEAR  OR 
C     QUADRATIC  CROSS  SECTION  FIT. 

SUBROUTINE  COEFF(KMAX.NORD) 

IMPLICIT  REAL*8  (A-H.O-Z) 


EXKERNEL.FOR  102 

06-25-1986 

COMMON/BLOCK1/  AA(3)  ,BB(3)  ,CC(3)  ,W(7)  ,SIG(7) 
C 
C     EVALUATE  COEFFICIENTS  FOR  LINEAR  FIT 

IF(NORD.EQ.l)  THEN 

DO  100  K=1,KMAX 

K1=K 

K2=K+1 
C 

AA(K)=0. 

BB(K)=(SIG(K1)-SIG(K2))/(W(K1)-W(K2)) 

CC(K)  =  (W(K1)*SIG(K2)-W(K2)*SIG(K1))/(W(K1)-W(K2)) 
100  CONTINUE 
C 
C     EVALUATE  COEFFICIENTS  FOR  QUADRATIC  FIT 

ELSE 

DO  200  K=1.KMAX 

K1=(K-1)*2+1 

K2=K1+1 

K3=Kl+2 
C 

D=W(K1)^2*(W(K2)-W(K3))-W(K2)^2*(W(K1)-W(K3))+W(K3)^2*(W(K1)-W( 
@K2)) 

D1^IG(K1)*(W(K2)-W(K3))-SIG(K2)*(W(K1)-W(K3))+SIG(K3)*(W(K1)-W(K2 
@)) 

D2=W(K1)^2*(SIG(K2)-SIG(K3))-W(K2)**2*(SIG(K1)-SIG(K3))+W(K3)**2* 
@(SIG(K1)-SIG(K2)) 

IB=W(K1)^2*(W(K2)*SIG(K3)-W(K3)*SIG(K2))-W(K2)**2*(W(K1)*SIG(K3)- 
@W(ra)*SIG(Kl))+W(K3)**2*(W(Kl)*SIG(K2)-W(K2)*SIG(Kl)) 
C 

AA(K)=D1/D 

BB(K)=D2/D 

CC(K)=D3/D 
200  CONTINUE 
C 

END  IF 
C 

500  RETURN 

END 
C 
C 
C 
C 

C     FUNCTION  SUBROUTINE  TO  EVALUATE  THE  EXACT  KERNEL  CROSS  SECTION 
C     BY  ANALYTICAL  INTEGRATION  OF  THE  PIECEWISE  CROSS  SECTION 
C     EXPANSION 

DOUBLE  PRECISION  FUNCTION  FNSIGG(I, J.KMAX) 

IMPLICIT  REAL*8(A-H.0-Z) 

DIMENSION  PHI (4) 

COMM0N/BL0CK1/  AA(3)  ,BB(3)  ,CC(3)  ,W(7)  ,SIG(7) 

COMMON/BL0CK2/  WB(4) ,SIGB(4) ,U(64) 

KMAX2=KMAX+1 

PI=4.*DATAN(1.D0) 

SIGMA=0. 
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Y=U(I)*U(J) 

Z=DSQRT(  ( 1 .  -U(  I)**2)*(  1 .  -U(  J)**2) ) 

DO  15  K=1.KMAX 

IF(U(I).EQ.-1..0R.U(I).EQ.1..0R.U(J).EQ.-1..0R.U(J).EQ.l.)   THEN 
DO  5  L=1,KMAX2 
IF(Y.EQ.WB(L))  THEN 
SIGMA=PI*SIGB(L) 
GO  TO  20 
END  IF 
5       CONTINUE 
K1=K+1 

DO  10  L=K,K1 
IF(WB(L).GT.Y)  THEN 

PHI(L)=0. 
ELSE 

PHI(L)=PI 
END  IF 
10   CONTINUE 
ELSE 

ARG1=DMIN1((WB(K)-Y)/Z,1.D0) 
ARG1 =DMAX1 ( - 1 . DO , ARG1 ) 
PHI(K)=DAC0S(ARG1) 
ARG2=DMIN1((WB(K+1)-Y)/Z,1.D0) 
ARG2=DMAX1(-1  .D0.ARG2) 
PHI(K+1)=DAC0S(ARG2) 
END  IF 

SIGMA^IGMA+((M(K)/40*Z^2)*(DSIN(2.*PHI(K))-DSIN(2.*PHI(K+1)))  + 
@(2.*M(K)*Y*Z+BB(K)*Z)*(DSIN(PHI(K))-DSIN(PHI{K+1)))  +  ((AA(K)/2.)*Z 
(^2MA(K)*Y**2+BB(K)*Y+CC(K))*(PHI(K)-PHI(K+1)) 
15  CONTINUE 

20  FNSIGG=2.*SIGMA 
RETURN 
END 
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C 

C  ***  PROGRAM  NAME:   ASKERNEL  *** 

C 

C     PROGRAM  TO  GENERATE  THE  APPROXIMATED  SCATTERING  KERNEL 

C     TRANSFER  MATRIX  FOR  A  GIVEN  DIFFERENTIAL  SCATTERING 

C     CROSS  SECTION.  DISCRETE-ORDINATES  SET,  AND  ORDER  OF 

C     ANGULAR  FLUX  EXPANSION.   THE  SCATTERING  CROSS  SECTION 

C     CAN  BE  REPRESENTED  BY  EITHER  A  PIECEWISE  LINEAR  OR 

C     PIECEWISE  QUADRATIC  EXPANSION  BETWEEN  THE  BREAKPOINTS 

C     IN  THE  SCATTERING  CROSS  SECTION.   THE  ANGULAR  FLUX  CAN  BE 

C     APPROXIMATED  BY  A  PIECEWISE  LINEAR  OR  PIECEWISE  QUADRATIC 

C     EXPANSION.   IN  ADDITION,  IF  A  PIECEWISE  LINEAR  EXPANSION  OF 

C     THE  SCATTERING  CROSS  SECTION  IS  USED,  THE  ANGULAR  FLUX  CAN 

C     BE  REPRESENTED  BY  A  PIECEWISE  CUBIC  EXPANSION. 

C 

C     INPUT  VARIABLES  ARE  AS  FOLLOWS: 

C     'NMAX-     =  ORDER  OF  FLUX  EXPANSION 

C     -M'        =  TOTAL  NUMBER  OF  FLUX  INTERVALS 

C     •NP'       =  TOTAL  NUMBER  OF  FLUX  DIRECTIONS  (NUMBER  OF 

C  DISCRETE  ORDINATES) 

C     •NBREAK•   =  NUMBER  OF  DISTINCT  BREAKPOINTS  IN  THE 

C  DIFFERENTIAL  SCATTERING  CROSS  SECTION 

C     'NORD'     =  ORDER  OF  CROSS  SECTION  EXPANSION 

C     •WB(I)'    =  BREAKPOINTS  IN  THE  SCATTERING  CROSS  SECTION 

C  (PLUS  THE  MIDPOINTS  OF  EACH  CROSS  SECTION 

C  SUBRANGE  IF  A  QUADRATIC  CROSS  SECTION 

C  EXPANSION  IS  USED) 

C     ■SIG(I)'   =  CROSS  SECTION  VALUES  AT  THE  WB(I)  VALUES 

C     'U(J)  '     =  DISCRETE  ORDINATES 

C 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  A(50,50)  ,SIG(7)  ,WB(7)  ,SIGT(50) 

COMMON/BLOCK1/  U(50)  ,W(4)  ,  AA(3)  ,BB(3) ,CC(3) 

C0MM0N/BL0CK2/  NMAX, KMAX.M.NP, NORD 
C 

C     READ  INPUT  DATA 
C 

READ(5,800)  NMAX.M.NP 
800  F0RMAT(3(I5)) 

READ(5,802)  NBREAK.NORD 

802  F0RMAT(2(I5)) 
C 

KMAX=NBREAK-1 
KMAX1=KMAX*N0RD+1 
NMAX1=NMAX+1 
LMAX=M+4 
C 

READ(5.803)  (WB(I)  ,  1=1  .KMAX1) 
READ(5,803)  (SIG(I)  ,  1=1  .KMAX1) 

803  F0RMAT(4(D18.8)) 
READ(5.804)  (U(J)  ,  J=l  ,NP) 

804  F0RMAT(3(D24.15)) 
C 

C     CALL  SUBROUTINE  TO  EVALUATE  ANGULAR  CROSS  SECTION  EXPANSION 
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C     COEFFICIENTS 
C 

CALL  XSEC(KMAX,NORD,SIG,WB) 
C 

DO  5  I=1,NBREAK 
W(I)=WB((I-l)*NORD+l) 
5  CONTINUE 
C 

C     CHECK  FOR  SYMMETRY  OF  THE  ANGULAR  MESH.   IF  THE  MESH 
C     IS  SYMMETRIC,  ONLY  ONE-HALF  OF  THE  APPROXIMATED 
C     SCATTERING  KERNEL  VALUES  A(I.J)  NEED  TO  BE  CALCULATED. 
C     THE  CITHERS  MAY  BE  DETERMINED  FROM  SYMMETRY  CONDITIONS. 
C 

IFLAG=1 
DO  10  1=1, NP 

IF(U(I).NE.-U(NP+1-I))  THEN 
IFLAG=0 
GO  TO  15 
END  IF 
10  CONTINUE 
15  CONTINUE 
C 

IF(IFLAG.EQ.l)  THEN 

NPl=IFIX((NP+l)/2.) 
ELSE 

NP1=NP 
END  IF 
C 

NP2=NP+1 
C 
C     BEGIN  LOOP  TO  CALCULATE  SCATTERING  KERNEL  MATRIX  ELEMENTS  A  (I,  J) 


C 


DO  20  1=1, NP1 
DO  20  J=1,NP 


A(I,J)=0.D0 
C 

C     CALL  SUBROUTINE  TO  DETERMINE  WHETHER  A  PARTICULAR  A(I.J) 
C     SCATTERING  KERNEL  IS  NONZERO. 
C 

CALL  C0VER(I,J,IFLAG2) 
C 

IF(IFLAG2.EQ.l)  THEN 
WRITE(9,920)  I.J 
920   F0RMAT(5X,I2.\  ' .  12) 
DO  25  N=1,NMAX1 
DO  25  L=1,LMAX 
DO  25  K=1.KMAX 
A(I,J)=A(I,J)+F(K,L,I,N,J) 
25   CONTINUE 
END  IF 
C 

A(I.J)=2.D0*A(I,J) 

IF(IFLAG.EQ.l)  A(NP2-I,NP2-J)=A(I, J) 
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20  CONTINUE 
C 

C     CHECK  ACCURACY  OF  A(I,J)  VALUES  BY  SUMMING  UP  THE  SCATTERING 
C     SOURCE  TERM  FOR  EACH  U(I),  ASSUMING  A  UNIFORM  FLUX.   IF  THE 
C     A(I,J)  VALUES  ARE  CORRECT.  THE  RESULTING  SOURCE  TERM  FOR 
C     EACH  U(I)  SHOULD  EQUAL  THE  GROUP  SCATTERING  CROSS  SECTION. 
C 

DO  50  1=1,  NP 

SIGT(I)=0. 

DO  50  J=1,NP 

SIGT(I)=SIGT(I)+A(J,I) 
50  CONTINUE 
C 

DO  60  1=1,  NP 

WRITE(6,900)  U(I),SIGT(I) 
60  CONTINUE 

900  F0RMAT(5X,  'SIGMAf  "  ,F9.6,  ' )  =  \D13.6) 
C 

WRITE(6,901) 

901  FORMAT(/////) 
C 

WRITE(6,902)  ((A(I, J) , J=l ,NP) , 1=1 .NP) 

902  F0RMAT(4(D18.8)) 
STOP 

END 
C 
C 
C 

C   XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXMXXXXMHXXXKXKHXXXXXXXXXXXXXXXXKXXXXX 

C  xxxxxxMXxxx»x»)()i)i)(x»)o<   SUBROUTINE  XSEC   xxxxxxxxioooooooooowxxxmcx 

C   XXXXXXXXXXXXXXXXXXXXXXXXXXXXXMXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C 

C     SUBROUTINE  TO  EVALUATE  EXPANSION  COEFFICIENTS  FOR  LINEAR  OR 

C     QUADRATIC  CROSS  SECTION  FIT. 

SUBROUTINE  XSEC(KMAX,NORD,SIG,WB) 
IMPLICIT  REAL*8  (A-H.O-Z) 
DIMENSION  SIG(7),WB(7) 

COMM0N/BL0CK1/  U(50)  ,W(4)  ,AA(3)  ,BB(3)  ,CC(3) 
C 

C     EVALUATE  COEFFICIENTS  FOR  LINEAR  FIT 
IF(NORD.EQ.l)  THEN 
DO  100  K=1,KMAX 
K1=K 
K2=K+1 
C 

AA(K)=0. 

BB(K)=(SIG(K1)-SIG(K2))/(WB(K1)-WB(K2)) 
CC(K)=(WB(K1)*SIG(K2)-WB(K2)*SIG(K1))/(WB(K1)-WB(K2)) 
100   CONTINUE 
C 

C       EVALUATE  COEFFICIENTS  FOR  QUADRATIC  FIT 
ELSE 

DO  200  K=1,KMAX 
K1=(K-1)*2+1 


ASKERNEL.FOR  107 

06-25-1986 

K2=K1+1 
K3=Kl+2 
C 

D=WB{Kl)x*2*(WB(K2)-WB(K3))-WB(K2)**2*(WB(Kl)-WB(K3))+WB(K3)**2* 
9  (WB(K1)-WB(K2)) 

Dl=SIG(Kl)*(WB(K2)-WB(ra))-SIG(K2)*(WB(Kl)-WB(K3))+SIG(K3)*(WB(K 
@  1)-WB(K2)) 

D2=WB(K1  )**2*(SIG(K2)-SIG(K3)  )-WB(K2)**2*(SIG(Kl  )-SIG(K3)  )+WB(K3 
0  )**2*(SIG(K1)-SIG(K2)) 

D3=WB(K1)**2*(WB{K2)*SIG(K3)-WB(K3)*SIG(K2))-WB(K2)**2*(WB(K1)*S 
@  IG(K3)-WB(K3)*SIG(K1))+WB(K3)**2*(WB(K1)*SIG(K2)-WB(K2)»SIG(K1)) 
C 

AA(K)=D1/D 
BB(K)=D2/D 
CC(K)=D3/D 
200   CONTINUE 
C 

END  IF 
C 

500  RETURN 
END 
C 
C 

c 

C       XXXXXXXXXXMXXXXXX>O<XXH)()000t)0(XXMXXX><)O(XXX><XXXXX><>OOOOOO<XXXXXXX>OOO<XXX 

C  xxxxxxxxxxxxxxxmxxmxxx   SUBROUTINE  COVER   xhxxxxxxxxxxxxxkmxxxxxxxx 

C       »><>OCXH)<XXXHXMXM)000000<MHXXXXXX)»CH)0000<XXKMX)00(X)C)00000<XXXH)0000<X)001KX 

C 

C     SUBROUTINE  TO  DETERMINE  WHETHER  A  PARTICULAR  SCATTERING  KERNEL 

C     A(I.J)  IS  NONZERO. 

SUBROUTINE  COVER(  I ,  J .  IFLAG2) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  THETA(50),UBAR(50) 

C0MM0N/BL0CK1/  U(50) ,W(4) ,AA(3) ,BB(3) ,CC(3) 

C0MM0N/BL0CK2/  NMAX.KMAX.M.NP.NORD 
C 

PI=4.*DATAN(1.D0) 

THEMAX=DAC0S(W(1)) 

THEMIN=DAC0S(W(KMAX+1 ) ) 
C 

DO  10  11=1, NP 

THETA(II)=DACOS(U(II)) 
10  CONTINUE 
C 

M1=M+1 

DO  20  JK=1,M1 

UBAR(JK)=U((JK-1)*NMAX+1) 
20  CONTINUE 
C 

IF(I.EQ.l)  THEN 
IL=1 
IU=NMAX+1 

ELSE  IF(I.EQ.NP)  THEN 
IL=NP-NMAX 
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IU=NP 
ELSE 

DO  30  JK=1,M1 
IF(U(I).EQ.UBAR(JK))  THEN 
IL=I-NMAX 
IU=I+NMAX 
GO  TO  40 
ELSE  IF(U(I).GT.UBAR(JK).AND.U(I).LT.UBAR(JK+1))  THEN 
IL=(JK-1)*NMAX+1 
IU=IL+NMAX 
GO  TO  40 
END  IF 
30   CONTINUE 
END  IF 
C 

C     DETERMINE  POLAR  BOUNDS  FOR  SCATTERING 
C 

40  CONTINUE 

DO  50  II=IL,IU 
IF(THETA(II)+THEMAX.LE.PI)  THEN 

BETMAX=THETA(  1 1  )+THEMAX 
ELSE  IF(THETA(II)+THEMIN.GE.PI)  THEN 

BETMAX=2 .  *PI-THETA(  1 1 )  -THEMIN 
ELSE 

BETMAX=PI 
END  IF 
C 

IF(THEMAX.LT.THETA(II))  THEN 

BETMIN=THETA(  1 1 ) -THEMAX 
ELSE  IF(THETA(II).LT. THEMIN)  THEN 

BETMIN=THEMIN-THETA( I I ) 
ELSE 

BETMIN=0. 
END  IF 
C 

IF(THETA(J).GE.BETMIN.AND.THETA(J).LE.BETMAX)  THEN 
IFLAG2=1 
GO  TO  500 
ELSE 

IFLAG2=0 
END  IF 
C 

50  CONTINUE 
C 

500  RETURN 
END 
C 
C 
C 

C  XXXMXXXXKXXXMXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXK 
C  XMXXXXMXXXXXXXXXXXXXXXXXXX  FUNCTION  F  xxxxxxxxxxxxxxxxxxxxxxxxxxx 
C   )00CXXXXXXXXXX)1XXXXXXX)0(XXXXXXXXXXXXXXX)C)CXXXXX)(XXXXXXXXXXXXKM)()0()1)0()(XX 

c 

C     SUBROUTINE  TO  EVALUATE  THE  FUNCTION  F(K,L,I,N,J) 
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DOUBLE  PRECISION  FUNCTION  F(K,L,  I.N.J) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON/BLOCK1/  U(50) ,W(4) ,AA(3) ,BB(3) ,CC(3) 

C0MM0N/BL0CK2/  NMAX.KMAX.M.NP.NORD 

COMMON/BLOCK3/  V(54) 
C 

C     CALL  SUBROUTINE  TO  EVALUATE  MU-PRIME  BREAKPOINTS  AND  ORDER  THE 
C     LIMITS  OF  INTEGRATION. 
C 

CALL  LIMITS(K.J) 
C 
C     CALCULATE  F(K,L, I.N.J) 

K1=K+1 

COEF=COEFF(L,N.I) 

IF(COEF.NE.O.DO)  THEN 
IF(NORD.EQ.l)  THEN 
C 

C       IF  A  LINEAR  CROSS  SECTION  EXPANSION  IS  USED,  THE  'SHORT' 
C       VERSION  OF  THE  FUNCTION  G(K,L,N,U,KK)  IS  UTILIZED. 
F=C0EF*(FNGA(K,L,N,U(J),K)-FNGA(K,L,N,U(J),K1)) 
ELSE 
C       IF  A  QUADRATIC  CROSS  SECTION  EXPANSION  IS  USED,  THE  'LONG' 
C       VERSION  OF  THE  FUNCTION  G(K,L,N.U,KK)  IS  UTILIZED 
F=eOEF*(FNGB(K,L,N,U(J),K)-FNGB(K.L,N.U(J),Kl)) 
END  IF 

ELSE 
F=O.DO 

END  IF 
C 

RETURN 

END 
C 
C 
C 

C   *><'"'''>'>»<»><  ><*"">»«XXXX>(XXXXXX)U0CX)t»»XXXH)0<XXX>O0<XXX)0(XXXX>O<XXXXmcxX)00O< 

C  *><>oooochxxxxxxx>o<xxxxx   SUBROUTINE  LIMITS   xxxxxxxxxxxxxxxiooootxxxx 

C   XXXXXX)H()<XXX)l)<XXX)()(XXX)()0(XXX)0()(XXXXXXXXXXXXXXX)0(XXXXX»X)(XX)<XXX)O()(X)()<X 

c 

C     SUBROUTINE  TO  EVALUATE  THE  MU-PRIME  BREAKPOINTS  AND  ORDER  THE 
C     INTEGRATION  LIMITS 

SUBROUTINE  LIMITS(K.J) 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  UPR(4),UBAR(50) 

COMMON/BLOCK1/  U(50) ,W(4) .AA(3) .BB(3) ,CC(3) 

C0MM0N/BL0CK2/  NMAX.KMAX.M.NP.NORD 

C0MM0N/BL0CK3/  V(54) 

EQUIVALENCE  (V(l)  ,UPR(1))  .  (V(5)  ,UBAR(1)) 
C 

M1=M+1 

M5=M+5 
C 

DO  10  JJ=1.M1 

UBAR(JJ)=U((JJ-1)*NMAX+1) 
10  CONTINUE 
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WW=W(K) 

WW1=W(K+1) 

UU=U(J) 

UPR(  1  )=UU*WW+DSQRT(  ( 1 .  -UU**2)*(  1 .  -WW**2)  ) 

UPR(2)=UU*WW-DSQRT( ( 1 . -UU**2)*( 1 . -WW**2) ) 

UPR(3)=UU*WW1+DSQRT(  ( 1 .  -UU**2)*{  1 .  -WW1**2) ) 

UPR(4)=UU»WW1-DSQRT(  ( 1 .  -UU**2)*(  1 .  -WW1**2)  ) 
C 

CALL  S0RT(V,M5) 
C 

RETURN 

END 
C 
C 

c 

C       XXXXXXXXXXXXXXXX>000<X><XXXXX»H»»XXXXXXXXXXXXXXXXXXXMHH)0(KXXXX»X)0000O<X 

C  xxxxxxxxxxxxxxxxxxxxxxxxx   SUBROUTINE  SORT   xxxxxxxxxxxxxxxxxxxxxxx 

C   XXXMX)(XXXXXXXXX)<X)()0()(XX>(XXXXX)<)<)0(XX)(XXX)(XXXXXX)(XXXX)(X)0()()(XXX)»0()0()CXXH 

c 

SUBROUTINE  SORT(A.N) 
IMPLICIT  REAL*8  (A-H.O-Z) 
DIMENSION  A(54) 
JUMP=N 
10  JUMP=JUMP/2 

IF(JUMP.EQ.O)  GO  TO  500 
J2=N-JUMP 
DO  25  J=1.J2 
15  I=J 
20  J3=I+JUMP 

IF(A(I).LE.A(J3))  GO  TO  25 
HOLD=A(I) 
A(I)=A(J3) 
A(J3)=HOLD 
I=I-JUMP 

IF(I.GT.O)  GO  TO  20 
25  CONTINUE 
GO  TO  10 
500  RETURN 
END 
C 
C 

c 

C   XX)()<XXXX>(XX)()()<XXXXXXKXXX)<>(XXXXXX)0()»(X)(XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C     Mxxxxxxxxxxxxxxxxioooooixx       FUNCTION  COEFF       xxxxxxxxxiohoodcxxkxxxxxx 

C       X)0(XXXXXXX)(X)()<XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)()()0(»)0(XXXXXXXXXXXXXX)»( 

c 

C     SUBROUTINE  TO  EVALUATE  THE  COEFFICIENTS  OF  THE  ANGULAR 
C     FLUX  EXPANSION. 

DOUBLE  PRECISION  FUNCTION  COEFF(L,N,I) 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  H(3).UBAR(50),UN(4) 

C0MM0N/BL0CK1/  U(50)  ,W(4)  ,AA(3)  ,BB(3)  ,CC(3) 

C0MM0N/BL0CK2/  NMAX.KMAX.M.NP.NORD 


AS KERNEL .FOR  m 

06-25-1986 


C0MM0N/BL0CK3/  V(54) 
C 

M1=M+1 

DO  100  JJ=1,M1 

UBAR(  JJ)=U(  (  JJ-1  )*NMAX+1 ) 
100  CONTINUE 
C 

DO  150  J=1,M 

H(1)=(U(I)-UBAR(J))*(UBAR(J+1)-U(I)) 

H(2)=(V(L)-UBAR(  J)  )*(UBAR(  J+l  )-V(L) ) 

H(3)=(V(L+1)-UBAR(J))*(UBAR(J+1)-V(L+1)) 


HH=1.D0 
C 

DO  200  11=1,3 
IF(H(II).GE.0.D0)  THEN 
H(II)=1.D0 
ELSE 

H(II)=O.D0 
END  IF 
C 

HH=HH*H(II) 
200  CONTINUE 
C 

O0EFF=0.D0 
C 

IF(HH.EQ.O.DO)  GO  TO  150 
C 

JL=J*NMAX-NMAX+1 
JU=JL+NMAX 
C 

IF(NMAX.EQ.2)  GO  TO  300 
IF(NMAX.EQ.3)  GO  TO  400 
C 

C     EVALUATE  EXPANSION  COEFFICIENTS  FOR  FIRST  ORDER  FIT 
C 

DO  250  K=JL.JU 
IF(U(K).NE.U(I))  THEN 
IF(N.EQ.l)  THEN 

COEFF=(-U(K)/(U(I)-U(K))) 
ELSE 

C0EFF=(1/(U(I)-U(K))) 
END  IF 
END  IF 
250  CONTINUE 

IF(OOEFF.NE.O.DO)  GO  TO  500 
GO  TO  150 
C 

300  CONTINUE 
C 

C     EVALUATE  EXPANSION  COEFFICIENTS  FOR  SECOND  ORDER  FIT 
C 

DEN0M=1.D0 
DNUM1=0.D0 


AS  KERNEL.  FOR  112 

06-25-1986 


DNUM2=1.D0 

DO  350  K=JL,JU 

IF(K.NE.I)  THEN 

DENOM=DENOM*(U(  I  )-U(K)  ) 
DNUM1=DNUM1-U(K) 
DNUM2=DNUM2*U(K) 
END  IF 
350  CONTINUE 
C 

IF(N.EQ.l)  THEN 

O0EFF=(DNUM2/DEN0M) 
ELSE  IF{N.EQ.2)  THEN 
OOEFF=(DNUM1/DENOM) 
ELSE 

COEFF=(1.DO/DENOM) 
END  IF 
C 

IF(COEFF.NE.O.DO)  GO  TO  500 
GO  TO  150 
C 

400  CONTINUE 
C 

C     EVALUATE  EXPANSION  COEFFICIENTS  FOR  THIRD  ORDER  FIT 
C 

DEN0M=1.D0 
DNUH1=0.D0 
DNUM2=0.D0 
DNUM3=1.D0 
JJ=0 

DO  450  K=JL,JU 
IF(K.NE.I)  THEN 
JJ=JJ+1 

DENOM=DENOM*(U( I )-U(K) ) 
DNUM1=DNUM1-U(K) 
DNUM3=DNUM3*U(K) 
UN(JJ)=U(K) 
END  IF 
450  CONTINUE 
C 

DNUM2=UN(  1  )*UN(2)+UN(  1  )*UN(3)+UN(2)*UN(3) 
C 

IF(N.EQ.l)  THEN 

C0EFF=-DNUM3/DEN0M 
ELSE  IF  (N.EQ.2)  THEN 

C0EFF=DNUM2/DEN0M 
ELSE  IF(N.EQ.3)  THEN 

C0EFF=DNUM1/DEN0M 
ELSE 

C0EFF=1.D0/DEN0M 
END  IF 
C 

IF(COEFF.NE.O.DO)  GO  TO  500 
C 

150  CONTINUE 
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C 

500  RETURN 
END 
C 

C 

c 

C       XX*MXXXXXX><XXXXK»XH)00<X)00<X)00<)000»0(HKXXHH>(XXXXX>00(XXXXX>OOOOOOOCHX)(XX 

C     mmxmmmxxhioooooocxmxxxxxxx       FUNCTION  FNGA       xxxxxxxxxxxxxxxxxxxxxxxxx 

C       ><XX)<X)<)<XXXKXX>OOOOO(MXXXH)OOO(X)0aO()0OOOCXXXXXXXXXMX><XXX)()»O<XXXXXXXXX)<XX 

c 

C     SUBROUTINE  TO  EVALUATE  THE  FUNCTION  FNG(K,L.N,U.KK) 
C     WHEN  A  PIECEWISE  LINEAR  CROSS  SECTION  EXPANSION  IS  USED 
DOUBLE  PRECISION  FUNCTION  FNGA(K.L,N,UU,KK) 
IMPLICIT  REAL*8  (A-H.O-Z) 

C0MM0N/BL0CK1/  U(50) ,W(4) ,AA(3) ,BB(3) ,CC(3) 
C0MH0N/BL0CK3/  V(54) 
C 

PI=4.D0*DATAN(1.D0) 
WW=W(KK) 
V1=V(L) 
V2=V(L+1) 
W=5.D-01*(V1+V2) 
C 

IF(WW.EQ.l.)  THEN 
FNGA=0. 
GO  TO  500 
END  IF 
C 

IF(V1.EQ.V2)  THEN 
FNGA=0. 
GO  TO  500 
END  IF 
C 

DNUM=WW-UU*W 

DEN0M=DSQRT(  ( 1 .  -UU**2)*(  1 .  -W**2)  ) 
C 

IF(DEN0M.EQ.0..0R.Vl.EQ.-1..0R.V2.EQ.l.)  THEN 
IF(DNUM.GT.O.)  THEN 
FNGA=0. 
GO  TO  500 
ELSE  IF(DNUM.LT.O.)  THEN 
ARG=-1. 
GO  TO  50 
ELSE  IF(DNUM.EQ.O.)  THEN 
IF(V1.EQ.-1.)  THEN 

Vl=-0. 99999 
END  IF 
IF(V2.EQ.1.)THEN 

V2=0. 99999 
END  IF 
GO  TO  100 
END  IF 
END  IF 
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40  CONTINUE 

ARG=DNUM/DENOM 

IF(ARG.LT.l.DO.AND.ARG.GT.-l.DO)  GO  TO  100 

50  CONTINUE 

IF(ARG.GE.l.DO)  THEN 
FNGA=O.DO 

ELSE  IF(ARG.LE.-l.DO)  THEN 

A1=(BB(K)*UU/(N+1))*(V2**(N+1)-V1**(N+1)) 

A2=(CC(K)/N)*(V2**N-V1**N) 

FNGA=PI*(A1+A2) 

END  IF 
GO  TO  500 

100  IF(N.EQ.l)  THEN 

FNGA=BB(K)*UU*(FNF1(V2,W,UU)-FNF1(V1,TO,UU))^E(K)*(FNF0(V2,WW, 
@  UU)-FNF0(V1,W,UU))+BB(K)*(FNG0(V2,W,UU)-FNG0(V1.WW,UU)) 

ELSE  IF(N.EQ.2)  THEN 

FNGA=BB(K)*UU*(FNF2(V2,W,UU)-FNF2(V1,TO,UU))+CC(K)*(FNF1(V2,WW, 
@  UU)-FNF1(V1,WW,UU))+BB(K)H(FNG1(V2,WW,UU)-FNG1(V1,WW,UU)) 

ELSE  IF(N.EQ.3)  THEN 

FNGA=BB(K)*UU*(FNF3(V2,WW,UU)-FNF3(V1,W.UU))+<X(K)*(FNF2(V2,WW, 
@  UU)-FNF2(V1,W,UU))+BB(K)*(FNG2(V2,WW,UU)-FNG2(V1,WW,UU)) 

ELSE  IF(N.EQ.4)  THEN 

FNGA=BB(K)»UUx(FNF4(V2,WW,UU)-FNF4(Vl,WW,UU))+CC(K)»(FNF3(V2,WW, 
@  UU)-FNF3(V1,W,UU))+BB(K)*(FNG3(V2,WW,UU)-FNG3(V1,WW,UU)) 

END  IF 


500  RETURN 
END 
C 

C 

c 

C   KM)0()()0(KX)()0<)()()<)l)0()()0()o<K)()0()()0O»<XXXXK)CX)»»0(K)()()0O()t)0()()0()()00000()0()»(XXX 

C  xx)(X)<)c)<)()(xx)i)<)<)()()oo(Xxxxx   FUNCTION  FNGB   xxxxxxxxxxxhkxhxxxmxxxxxxx 

C   X)<X)()»0(X)C)(»)()()()<X)()()0()0»<XX)(»)(X)(XXXX)()()0()(XXXKX)<X)()(XXX)»()()()(XX)»()(HXX)(XXXX 

c 

C     SUBROUTINE  TO  EVALUATE  THE  FUNCTION  FNG(K.L,N,U,KK) 

C     WHEN  A  PIECEWISE  QUADRATIC  CROSS  SECTION  EXPANSION  IS  USED 

DOUBLE  PRECISION  FUNCTION  FNGB(K,L.N,UU,KK) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON/BLOCKl/  U(50)  ,W(4)  ,AA(3)  ,BB(3)  ,CC(3) 

C0MM0N/BL0CK3/  V(54) 
C 

PI=4.D0*DATAN(1.D0) 

WW=W(KK) 
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V1=V(L) 

V2=V(L+1) 

W=5.D-01*(V1+V2) 

IF(WW.EQ.l.)   THEN 

FNGB=0. 

GO  TO  500 
END  IF 

IF(V1.EQ.V2)  THEN 

FNGS=0. 

GO  TO  500 
END  IF 

DNUM=WW-UU*W 

DENOM=DSQRT(  ( 1 .  -UU**2)*(  1 .  -W**2)  ) 

IF(DENOM.EQ.0..OR.Vl.EQ.-l..OR.V2.EQ.l.)  THEN 
IF(DNUM.GT.O.)  THEN 
FNGB=0. 
GO  TO  500 
ELSE  IF(DNUM.LT.O.)  THEN 
ARG=-1. 
GO  TO  50 
ELSE   IF(DNUM.EQ.O.)  THEN 
IF(V1.EQ.-1.)  THEN 

Vl=-0. 99999 
END  IF 
IF(V2.EQ.l.)  THEN 

V2=0. 99999 
END  IF 
END  IF 
END  IF 

40  CONTINUE 

ARG=DNUM/DEN0M 

IF(ARG.LT.l.DO.AND.ARG.GT.-l.DO)  GOTO  100 

50  CONTINUE 

IF(ARG.GE.l.DO)  THEN 
FNGB=0.DO 

ELSE  IF(ARG.LE.-l.DO)  THEN 

Al=(AA(K)*UU**2/(N+2)-AA(K)*( 1 . -UU**2)/(2 . *(N+2) ) )*(V2**(N+2)-Vl 
@  **(N+2)) 

A2=(BB(K)*UU/(N+1 )  )*(V2**(N+1  )-yi**(N+l)  ) 
A3=(CC(K)/N+AA(K)*(  1 .  -UU**2)/(2-  *N)  )*(V2**N-V1**N) 
FNGB=PI*(A1+A2+A3) 

END  IF 
GO  TO  500 
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100  IF(N.EQ.l)  THEN 

FNGB=M(K)/2.*(3.*UU**2-10K™F2(V2,WW,UU)-FNF2(V1,WW,UU))+BB(K 
@  )*UU«(FOTl(V2,WW,UU)-FNFl(Vl,W,UU))  +  ((r(K)+M(K)*(l.-UU**2)/2.) 
@  »«(FNF0(V2,WW.UU)-FNF0(Vl,WW,UU))+3./2.»*AA(K)«UU*«(FNGl(V2,WW,UU)- 
@     FNG1(V1.WW,UU))+(BB(K)+AA(K)«WW/2.)*«(FNG0(V2,WW,UU)-FNG0(V1.WW,U 

9    u)) 

c 

ELSE  IF(N.EQ.2)  THEN 

FNGB=M(K)/2.*(3.*UU**2-10*(FNF3(V2,WW,UU)-FNF3(V1,WW,UU))+BB(K 
@  )*UU*^FNF2(V2,W,UU)-FNF2(Vl,W,UU))  +  (CC(K)+M(K)*(l.-UU**2)/2.) 
@  *(FNF1(V2,OT,UU)-FNF1(V1,W,UU))+3V2.*AA(K)*UU*(FNG2(V2,WW,UU)- 
0     FTJG2(Vl,W,UU))+(BB(K)+M(K)xWW/2.)*(FNGl(V2,WW,UU)-FNGl(Vl,WW,U 

e  u)) 

c 

ELSE  IF(N.EQ.3)  THEN 

FNGB=AA(K)/2.»(3.s«UU^2-l.)'»(FNF4(V2,WW.UU)-FNF4(Vl,WW,UU))+BB(K 
@  )*UU*(FNF3(V2.TO.UU)-FNF3(Vl,W,UU))+(CC(K)+M(K)*(l.-UU**2)/2.) 
@  *(FNF2(V2,W,UU)-FNF2(Vl,WW,UU))+3./2.xAA(K)xUU*(FNG3(V2,WW.UU)- 
@  FNG3(V1.W,UU))  +  (BB(K)+M(K)*WW/2.)*(FNG2(V2,WW,UU)-FNG2(V1,WW,U 
@  U» 
END  IF 
C 

500  RETURN 
END 
C 
C 
C 
C 

C   XMXXXXXXXXXXXXXXXKXIOOCKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C  xxii)i)[)()(i[)[)[K)i oodODOi   FUNCTION  SUBROUTINES   xxxxxxxxxxxxxxxxxxxxxx 

C   XXXXXXXXXHHXXXXXXXXXXXXXXXXXXXXKXXHKKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

c 

C     SUBROUTINE  TO  EVALUATE  THE  FUNCTION  F0(X,W,U) 

DOUBLE  PRECISION  FUNCTION  FNFO(X.W.U) 

IMPLICIT  REAL*8  (A-H.O-Z) 

A=W/DSQRT(1.-U**2) 

B=U/DSQRT(1.-U**2) 

A1=A/DSQRT(B**2+1 .  ) 

FNF0=X*FNH1(X,A,B)+A1*FNH2(X,A,B)+FNH3(X,A,B)-FNH4(X,A,B) 

RETURN 

END 
C 
C 
C 
C 
C     SUBROUTINE  TO  EVALUATE  THE  FUNCTION  F1(X,W,U) 

DOUBLE  PRECISION  FUNCTION  FNF1(X,W,U) 

IMPLICIT  REAL*8  (A-H.O-Z) 

A=W/DSQRT(1.-U**2) 

B=U/DSQRT(  1 .  -U**2) 

A1=(A**2*B-B*(B**2+1 .  ) )/(  (DSQRT(B**2+1 .  )  )**3) 

FI^1=(X^2*FNH1(X,A,B)+A1*FNH2(X,A.B)-FNH3(X,A,B)-FNH4(X,A.B)+A*FN 
@H5(X,A,B))/2. 

RETURN 
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END 
C 

c 
c 
c 

C     SUBROUTINE  TO  EVALUATE  THE  FUNCTION  F2(X,W,U) 

DOUBLE  PRECISION  FUNCTION  FNF2(X,W,U) 

IMPLICIT  REAL*8  (A-H.O-Z) 

A=W/DSQRT(1.-U**2) 

B=U/DSQRT(1.-U**2) 

Al=3.*A**2*B**2/(  (B**2+l .  )**2) 

A2=Al+2.-(2.*B**2+A**2-l .  )/(B**2+l .  ) 

A3=A*A2/(2.*DSQRT(B**2+1.)) 

A4=  (  A*X-2 . *B+ ( 3 .  *A**2*B/  {  B**2+ 1 .  )  )  )  /2 . 

FNF2=(X^3*FNH1(X.A,B)M3*FNH2(X,A,B)+FNH3(X,A,B)-FNH4(X,A,B)+A4*F 
@NH5(X,A,B))/3. 

RETURN 

END 
C 
C 
C 
C 
C     FUNCTION  SUBROUTINE  TO  EVALUATE  F3(X,W,U) 

DOUBLE  PRECISION  FUNCTION  FNF3(X,W.U) 

IMPLICIT  REAL*8( A-H.O-Z) 

A=W/DSQRT(1.-U**2) 

B=U/DSQRT(1.-U**2) 

P=1./(B**2+1.) 

R=DSQRT(B**2+1 . ) 

Al=(5 .  *A**2*B**3*P-3 .  *B*(  A**2-l .  )  )*A**2/(2 .  *R**5) 

A2=  (  3 .  *A**2*B**2*P-3 .  *A**2+ 1 .  +2 .  *R**2  )  *B/  (  2 .  *R**3  ) 

A3=(X+3 . *A*B*P)*B/2 . 

A4=(X**2/3 . +5 . *A*B*P*X/6 . -2 . *(A**2-1 . )*P/3 . +1 . ) 

A5=A*(  A4+5 .  *A**2*B*«2*P**2/2  .  ) 

FOT3=(X^*™i(X,A,B)  +  (Al-A2)*FNH2(X,A,B)-FNH3(X,A,B)-FNH4(X,A,B) 
@+(A5-A3)*FNH5(X,A,B))/4. 

RETURN 

END 
C 
C 
C 
C 
C     SUBROUTINE  TO  EVALUATE  THE  FUNCTION  F4(X,W,U) 

DOUBLE  PRECISION  FUNCTION  FNF4(X,W,U) 

IMPLICIT  REAL*8  (A-H.O-Z) 

A=W/DSQRT( 1 . -U**2) 

B=U/DSQRT(1.-U**2) 

P=1./(B**2+1.) 

R=DSQRT(B**2+1.) 

S=A**2-1. 

T=A*B 

A 1 = (  7 .  *A*T/  (  4 .  *R**2  )  -B  )  *  (  5 .  *T**2*P-3 .  *S )  *T/  (  2 .  *R**5  ) 

A2=(3 .  *A*T**2*P/2 .  -A*S/2.  -B*T+A*R**2)/(R**3) 

A3=(3 .  *A*S/(8 .  *R**5)  )*(3 .  *T**2*P-S) 
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A2=A1+A2-A3 

A5= (X**2/ ( 3 . *P } +5 . *T*X/6 . -2 . *S/3 . +5 . *T**2* P/2 . ) *7 . *A*T*P**2/4 . 

A6= ( X+3 . *T*P+ ( X**3 ) /2 . ) *A/2 . - ( X+3 .  *T*P ) *3 . *A*S*P/8 . 

A7= (  X**2/3 .  +5 .  *T*P*X/6  .  -2 .  *S*P/3 .  +5 .  *T**2*P**2/2 . + 1 . )  *B 

A5=A5+A6-A7 

FNF4=(X^5*FNH1(X,A,B)+A2*FNH2(X,A,B)+FNH3(X.A.B)-FNH4(X,A,B)+A5*F 
@NH5(X.A.B))/5. 

RETURN 

END 
C 
C 
C 
C     SUBROUTINE  TO  EVALUATE  THE  FUNCTION  GO(X,W,U) 

DOUBLE  PRECISION  FUNCTION  FNGO(X.W.U) 

IMPLICIT  REAL*8  (A-H.O-Z) 

C=l.-U**2-W**2 

D=2.*W*U 

ARG1=C+D»X-X**2 

IF(ARGl.LT.0.D0.AND.ARGl.GT.-l.D-05)  ARG1=0.D0 

Al=(2 .  *X-D)*DSQRT(  ARG1  )/4 . 

A2=(D**2+4.*C)/8. 

FNG0=A1-A2*FNH6(X,C,D) 

RETURN 

END 
C 
C 
C 
C 
C     SUBROUTINE  TO  EVALUATE  THE  FUNCTION  Gl(X.W.U) 

DOUBLE  PRECISION  FUNCTION  FNGl(X.W.U) 

IMPLICIT  REAL*8  (A-H.O-Z) 

C=l.-U**2-W**2 

D=2.*W*U 

ARG1=C+D*X-X**2 

IF(ARG1.LT.O.DO.AND.ARG1.GT.-1.D-05)  ARGl=O.DO 

Al=( (DSQRT(ARG1 ) )**3)/3. 

A2=DSQRT(  ARG1  )*D*(D-2 .  *X)/8 . 

A3=D*(D**2+4.*C)/16. 

FNG1=-A1-A2-A3*FNH6(X,C,D) 

RETURN 

END 
C 
C 
C 
C 
C     SUBROUTINE  TO  EVALUATE  THE  FUNCTION  G2(X.W,U) 

DOUBLE  PRECISION  FUNCTION  FNG2(X,W,U) 

IMPLICIT  REAL*8  (A-H.O-Z) 

C=l.-U**2-W**2 

D=2.*W*U 

ARG1=C+D*X-X**2 

IF(ARGl.LT.0.D0.AND.ARGl.GT.-l.D-05)  ARGl=O.DO 

Al=(  (DSQRT(ARG1 )  )**3)/4. 

A2=Al*(X+5.*D/6.) 
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A3=(5 . *D»*2+4 .  *C)/16 . 
FNG2=-A2+A3*FNG0(X , W , U) 
RETURN 
END 
C 

c 
c 

c 

C  SUBROUTINE  TO  EVALUATE  THE  FUNCTION  G3(X,W,U) 

DOUBLE  PRECISION  FUNCTION  FNG3(X,W.U) 

IMPLICIT  REAL*8   (A-H.O-Z) 

C=l .  -U**2-W**2 

D=2.*W*U 

ARG1=C+D*X-X**2 

IF(ARGl.LT.0.D0.AND.ARGl.GT.-l.D-05)   ARG1=0.D0 

Al={  (DSQRT(ARG1 )  )**3)*X**2/5 . 

A2=7.*D/10. 

A3=2.«C/5. 

FNG3=-A1+A2*FNG2(X,W,U)+A3*FNG1(X,W,U) 

RETURN 

END 
C 
C 
C 
C 
C  FUNCTION  SUBROUTINE  FOR  Hl(X.A.B) 

DOUBLE  PRECISION  FUNCTION  FNH1(X,A,B) 

IMPLICIT  REAL*8   (A-H.O-Z) 

ARG1=(A-B*X)/DSQRT(  1 .  -X**2) 

IF(ARG1 .GT. 1 .DO. AND.DABS(ARG1-1 .DO) .LT. 1 .D-05)  ARG1=1 .DO 

IF(ARGl.LT.-l.D0.AND.DABS(ARGl+l.D0).LT.l.D-05)  ARG1=-1.D0 

FNH1=DARC0S(ARG1) 

RETURN 

END 
C 
C 
C 
C 
C     FUNCTION  SUBROUTINE  FOR  H2(X,A,B) 

DOUBLE  PRECISION  FUNCTION  FNH2(X.A,B) 

IMPLICIT  REAL*8  (A-H.O-Z) 

Al=-1 .  D0*(B**2+1  .  DO)*X+A*B 

A2=DSQRT( 1 . D0-A**2+B**2) 

ARG1=A1/A2 

IF(ARGl.GT.l.D0.AND.DABS(ARGl-l.D0).LT.l.D-05)  ARG1=1.D0 

IF(ARGl.LT.-l.D0.AND.DABS(ARGl+l.D0).LT.l.D-05)  ARG1=-1.D0 

FNH2=DARSIN(ARG1) 

RETURN 

END 
C 
C 

c 
c 

C     FUNCTION  SUBROUTINE  TO  EVALUATE  H3(X,A,B) 
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DOUBLE  PRECISION  FUNCTION  FNH3(X,A,B) 

IMPLICIT  REAL*8  (A-H.O-Z) 

IF(A.EQ.-B)  THEN 
FNH3=0. 
GO  TO  500 

END  IF 

A1=(A*B+B**2+1.D0)*(X+1.D0)-((A+B)**2) 

A2=DABS  ( X+ 1 )  *DSQRT  ( 1 .  D0-A**2+B**2 ) 

A3=5 .  D-01*(A+B)/DABS(  A+B) 

ARG1=A1/A2 

IF(ARG1.GT.1.DO.AND.DABS(ARG1-1.DO).LT.1.D-05)  ARG1=1.D0 

IF(ARGl.LT.-l.D0.AND.DABS(ARGl+l.D0).LT.l.D-05)  ARG1=-1.D0 

FNH3=A3*DARSIN(ARG1 ) 
500  RETURN 

END 
C 
C 
C 
C 
C     FUNCTION  SUBROUTINE  TO  EVALUATE  THE  FUNCTION  H4(X.A,B) 

DOUBLE  PRECISION  FUNCTION  FNH4(X,A,B) 

IMPLICIT  REAL*8  (A-H.O-Z) 

IF(A.EQ.B)  THEN 
FNH4=0. 
GO  TO  500 

END  IF 

A1=(A*B-B**2-1 .D0)*(X-1 .D0)-( (A-B)**2) 

A2=DABS  (  X- 1 .  DO )  *DSQRT  ( 1 .  D0-A**2+B**2  ) 

A3=5 .  D-01*(A-B)/DABS(A-B) 

ARG1=A1/A2 

IF(ARG1 .GT. 1 .DO. AND.DABS(ARG1-1 .DO) .LT. 1 .D-05)  ARG1=1 .DO 

IF(ARGl.LT.-l.D0.AND.DABS(ARGl+l.D0).LT.l.D-05)  ARG1=-1.D0 

FNH4=A3*DARSIN(  ARG1 ) 
500  RETURN 

END 
C 
C 

c 

C 

C     FUNCTION  SUBROUTINE  TO  EVALUATE  THE  FUNCTION  H5(X,A,B) 

DOUBLE  PRECISION  FUNCTION  FNH5(X,A,B) 

IMPLICIT  REAL*8  (A-H.O-Z) 

ARG1=(2.D0*A*B*X-(B**2+1  .D0)*X**2-(A**2-1  .DO)) 

IF(ARG1.LT.O.DO.AND.ARG1.CT.-1.D-05)  ARG1=0.D0 

A1=DSQRT(ARG1) 

A2=B**2+1.D0 

FNH5=A1/A2 

RETURN 

END 
C 
C 
C 
C 
C     SUBROUTINE  TO  EVALUATE  THE  FUNCTION  H6(X.C,D) 
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DOUBLE  PRECISION  FUNCTION  FNH6(X,C,D) 

IMPLICIT  REAL*8  (A-H.O-Z) 

ARGl=(D-2 .  *X)/DSQRT(D**2+4 .  *C) 

IF(ARG1  .CT.  1  .DO. AND.DABS(ARG1-1  .DO)  .LT.  1  .D-05)  AKG1-1  .DO 

IF(ARGl.LT.-l.D0.AND.DABS(ARGl+l.D0).LT.l.D-05)  ARG1=-1.D0 

FNH6=DARSIN(ARG1) 

RETURN 

END 
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C 

C  ***     PROGRAM  NAME:      TRANS     *** 

C 

C  PROGRAM  TO  SOLVE  THE  ONE-DIMENSIONAL.  AZIMUTHALLY  SYMMETRIC 

C  PLANE  GEOMETRY  DISCRETE  ORDINATES  EQUATIONS.   SCATTERING 

C  CAN  BE  ISOTROPIC  OR  ANISOTROPIC.   ANISOTROPIC  GROUP-TO-GROUP 

C  APPROXIMATED  SCATTERING  KERNELS  ARE  GENERATED  BY  THE  CODE 

C  'ASKERNEL'  AND  THE  RESULTING  APPROXIMATED  SCATTERING  KERNEL 

C  TRANSFER  MATRIX  IS  INPUT  INTO  THIS  CODE.   BECAUSE  THIS  CODE 

C  USES  THE  APPROXIMATED  SCATTERING  KERNEL  TECHNIQUE  FOR 

C  EVALUATION  OF  THE  SCATTERING  SOURCE  TERM,  THE  DISCRETE 

C  ORDINATES  CAN  BE  ARBITRARILY  CHOSEN  (I.E.  ,  THEY  DO  NOT  HAVE 

C  TO  BE  A  STANDARD  NUMERICAL  QUADRATURE  SET) . 

C 

C 

C  MAJOR  VARIABLES  ARE  AS  FOLLOWS: 

C  'NR'      =  NUMBER  OF  REFLECTED  DISCRETE  DIRECTIONS 

C  'NT'      =  NUMBER  OF  TRANSMITTED  DISCRETE  DIRECTIONS 

C  •NTCIT•    =  TOTAL  NUMBER  OF  DISCRETE  ORDINATES 

C  'M'       =  NUMBER  OF  SPATIAL  NODES.   THE  NODES  ARE  DEFINED 

C  TO  BE  AT  THE  EDGES  OF  THE  MESH  CELLS. 

C  •WIDTH'   =  WIDTH  OF  SLAB 

C  'DEL'     =  WIDTH  OF  EACH  SPATIAL  MESH  CELL 

C  'IFLAG1'   =  BOUNDARY  CONDITION  FLAG.   SET  IFLAG1=0  FOR  A  UNIT 

C  FLUX  BOUNDARY  CONDITION.   SET  IFLAG1=1  FOR  A  UNIT 

C  CURRENT  BOUNDARY  CONDITION. 

C  'IFLAG2'   =  UNITS  FLAG.   SET  IFLAG2=0  FOR  UNITS  OF  CENTIMETERS. 

C  SET  IFLAG2=1  FOR  UNITS  OF  MEAN  FREE  PATHS. 

C  'NEXP'    =  ORDER  OF  FLUX  EXPANSION 

C  •NINT'    =  NUMBER  OF  FLUX  INTERVALS 

C  'EPS'     =  CONVERGENCE  CRITERIA  FOR  THE  ANGULAR  FLUXES 

C  'IS'      =  SOURCE  ANGLE  (DISCRETE  DIRECTION  NUMBER) 

C  'SIGT'    =  TOTAL  GROUP  CROSS  SECTION 

C  'CAP'     =  CONVERGENCE  ACCELERATION  PARAMETER 

C  'NMAX'    =  MAXIMUM  NUMBER  OF  ITERATIONS  TO  BE  PERFORMED 

C  'NUMB'    =  NUMBER  DENSITY  OF  SCATTERING  MATERIAL 

C  'SIGGT'   =  TOTAL  GROUP  SCATTERING  CROSS  SECTION  (NOT  TO  BE 

C  CONFUSED  WITH  THE  TOTAL  GROUP  CROSS  SECTION) 

C 

C  NOTE:  THE  ORDER  OF  THE  FLUX  EXPANSION  ('NEXP')  CAN  BE 

C  1  (PIECEWISE  LINEAR),  2  (PIECEWISE  QUADRATIC),  OR 

C         3  (PIECEWISE  CUBIC) .  THE  TOTAL  NUMBER  OF  DISCRETE 

C  ORDINATES  USED  MUST  EQUAL  (NEXP*NINT)+1 .   THE  VALUES 

C  OF  THESE  VARIABLES  MUST  BE  THE  SAME  AS  THE  VALUES 

C  USED  IN  THE  CODE  •ASKERNEL'  TO  GENERATE  THE 

C  APPROXIMATED  SCATTERING  KERNEL  TRANSFER  MATRIX. 

C 

C  NOTE:  THE  INPUT  VALUES  OF  'SIGT',  'NUMB',  AND  'SIGGT'  DEPEND  ON 

C  WHETHER  UNITS  OF  CENTIMETERS  OR  MEAN  FREE  PATHS  ARE  USED. 

C  THESE  VALUES  SHOULD  BE  INPUT  AS  FOLLOWS: 

C 

C  CASE  (1):  IF  UNITS  OF  CENTIMETERS  ARE  USED  (IFLAG2=0), 

C  INPUT  THE  TOTAL  GROUP  MICROSCOPIC  CROSS  SECTION 

C  (IN  BARNS)  FOR  'SIGT'  (THIS  INCLUDES  SCATTERING 
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C  AND  ALL  OTHER  INTERACTIONS).   INPUT  THE  NUMBER 

C  DENSITY  (1/CM**3)  OF  THE  SCATTERING  MATERIAL  FOR 

C  •NUMB',  AND  INPUT  THE  TOTAL  GROUP  MICROSCOPIC 

C  CROSS  SECTION  (IN  BARNS)  FOR  ,SIGGT'. 

C 

C  CASE  (2):  IF  UNITS  OF  MEAN  FREE  PATHS  ARE  USED  (IFLAG2=1), 

C  INPUT  A  VALUE  OF  '  1 '  FOR  •SIGT'  AND  THE  VALUE  OF 

C  C  (THE  MEAN  NUMBER  OF  SECONDARIES  PER  COLLISION) 

C  FOR  'NUMB'.   INPUT  THE  TOTAL  GROUP  SCATTERING 

C  CROSS  SECTION  (IN  BARNS)  FOR  •SIGGT1  . 

C 

C  THE  FOLLOWING  ARRAYS  ARE  UTILIZED: 

C  F(M.N)   =  NEUTRON  FLUX  AT  CELL  EDGE 

C  G(M.N)   =  NEUTRON  FLUX  AT  CELL  MIDPOINT 

C  Q(M,N)   =  INSCATTER  SOURCE 

C  A(I.J)   =  APPROXIMATED  SCATTERING  KERNEL  TRANSFER  MATRIX 

C  U(N)     =  ANGULAR  MESH  POINTS  (THE  DISCRETE  ORDINATES) 

C 

IMPLICIT  REAL*8(A-H,0-Z) 

REAL*8  NORM, NUMB 

COMMON/B1/  F(200.50),G(200,50),Q(200,50),U(50),A(50,50) 

C0MM0N/B2/  M.NR,NT,IS,WIDTH,ERR1,ERR2,ERR3 
C 

C     READ  INPUT  PARAMETERS  (DEFINED  ABOVE) 
C 

READ(5,800)  M, WIDTH,  IFLAG1 ,  IFLAG2 

800  F0RMAT(I5,D20.4,5X,2(I5)) 
READ(5.803)  NEXP.NINT 
READ(5,801)  EPS 

801  F0RMAT(5X,D20.4) 
READ(5.802)  IS.NMAX.SIGT.CAP 

802  F0RMAT(2(I5),D17.6,D18.4) 
READ(5,803)  NR.NT 

803  F0RMAT(2(I5)) 
C 

MM=M-1 

XMM=DFLOAT(MM) 
DEL=WIDTH/XMM 
NTOT=NR+NT 
C 

READ(5.804)  (U(I)  ,  1=1  ,NTOT) 

804  FORMAT(3(D24.10)) 
C 

READ(5.805)  NUMB.SIGGT 

805  F0RMAT(2(D24.10)) 
C 

IF(IFLAG2.EQ.O)  THEN 
CC=SIGGT/SIGT 
SIGT=SIGT*NUMB 
SIGGT2=1 . 
ELSE 
CC=NUMB 
SIGGT2=SIGGT 
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END  IF 
C 
C  READ  IN  APPROXIMATED  SCATTERING  KERNEL  TRANSFER  MATRIX 

READ(5.806)    ((A(I, J) , J=l ,NTOT) , 1=1 ,NTOT) 
806  F0RMAT(4(D18.8)) 
C 

DO  10  I=1,NT0T 

DO  10  J=1,NT0T 

A(  I ,  J)=A(  I ,  J)*NUMB/SIGGT2 
10  CONTINUE 
C 

C     INITIALIZE  FLUX  AND  SOURCE  ARRAYS  AT  ZERO 
C 

DO  40  1=1,  M 

DO  45  J=1,NT0T 

F(I.J)=O.DO 

Q(I,J)=0.D0 
45  CONTINUE 
40  CONTINUE 
C 

C     SET  INCIDENT  FLUX  USING  BOUNDARY  CONDITION 
C 

CALL  N0RMAL(U,IS,IFLAG1,N0RM,NEXP,NINT) 
C 

F(1,IS)=1./N0RM 
C 

C     CALL  SUBROUTINE  TO  SOLVE  D.O.  EQUATIONS 
C 

CALL  DO(DEL, EPS, SIGT, CAP, NMAX,  ITER) 
C 
C 

C     CALL  SUBROUTINE  TO  PRINT  OUT  RESULTS 
C 


CALL  0UTPUT(IFLAG1,IFLAG2,  ITER, NMAX, EPS, SIGT, CC,NEXP,NINT) 


STOP 

END 
C 
C 
C 
C     SUBROUTINE  TO  CALCULATE  SOURCE  NORMALIZATION 

SUBROUTINE  NORMAL(U,  II.IFLAGl.NORM.NEXP.NINT) 

IMPLICIT  REAL*8  (A-H.O-Z) 

REAL*8  NORM 

DIMENSION  U(50) 

NORM=0. 

DO  50  J=1,NINT 

N0RM=N0RM+S0RCE(  J,  II  ,U,  IFLAG1  .NEXP.NINT) 
50  CONTINUE 

RETURN 

END 
C 
C 
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DOUBLE  PRECISION  FUNCTION  SORCEfJ,  II, U,  IFLAG1  ,NEXP,NINT) 
IMPLICIT  REAL*8   (A-H.O-Z) 
DIMENSION  U(50),UBAR(50),UN(3) 
NINT1=NINT+1 
DO  10  JJ=1,NINT1 
UBAR(JJ)=U((JJ-1)*NEXP+1) 
10  CONTINUE 

H=(U(II)-UBAR(J))*(UBAR(J+1)-U(II)) 
SORCE=0. 

IF(H.LT.O.)  GO  TO  500 
C 

JL=J*NEXP-NEXP+1 
JU=JL+NEXP 
UJ1=UBAR(J+1) 
UJ=UBAR(J) 

IF(NEXP.EQ.2)  GO  TO  450 
IF(NEXP.EQ.3)  GO  TO  480 
C 

C     EVALUATE  NORMALIZATION  FOR  LINEAR  FLUX  EXPANSION 
•  DO  400  K=JL,JU 
IF(K.NE.II)  THEN 

IF(IFLAGl.EQ.O)  THEN 

SORCE=(  (UJl**2-UJ**2)/2.  -(UJ1-UJ)*U(K)  )/(U(  II  )-U(K)  ) 
ELSE  IF(IFLAGl.EQ.l)  THEN 

S0RCE=(  (UJl**3-UJ**3)/3 .  -(UJl**2-UJ**2)*U(K)/2 .  )/(U(  1 1  )-U(K)  ) 
END  IF 
GO  TO  500 
END  IF 
400  CONTINUE 
C 

450  CONTINUE 
C     EVALUATE  NORMALIZATION  FOR  QUADRATIC  FLUX  EXPANSION 
DEN0M=1.D0 
DN1=0.D0 
DN2=1.D0 
DO  475  K=JL,JU 
IF(K.NE.II)  THEN 

DEN0M=DEN0M*(U(  II)-U(K)  ) 
DN1=DN1+U(K) 
DN2=DN2*U(K) 
END  IF 
475  CONTINUE 

IF(IFLAGl.EQ.O)  THEN 

S0RC^((UJl**3-UJ*w3)/3.-(UJl**2-UJ**2)*DNl/2.+(UJl-UJ)*DN2)/DEN 
@  OM 
ELSE  IF(IFLAGl.EQ.l)  THEN 

S0RC£=((UJ1^4-UJ>^)/4.-(UJl>«3-UJ**3)*DNl/3.  +  (UJl**2-UJ**2)*DN 
@  2/2. )/DENOM 
END  IF 
GO  TO  500 
C 

480  CONTINUE 
C     EVALUATE  NORMALIZATION  FOR  CUBIC  FLUX  EXPANSION 
DEN0M=1.D0 
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DN1=0.D0 
DN2=0.D0 
DN3=1.D0 
JJ=0 

DO  485  K=JL,JU 
IF(K.NE.II)  THEN 
JJ=JJ+1 

DEN0M=DEN0M*(U(  1 1  )-U(K)  ) 
DN1=DN1+U(K) 
DN3=DN3*U(K) 
UN(JJ)=U(K) 
END  IF 
485  CONTINUE 

DN2=UN(1)*UN(2)+UN(1)*UN(3)+UN(2)*UN(3) 
IF(IFLAGl.EQ.O)  THEN 

SORCE=(  (UJl**4-UJ**4)/4.  -(UJl**3-UJ**3)*DNl/3 .  +  (UJ1**2-UJ**2)*DN 
9     2/2.-(UJl-UJ)*DN3)/DEN0M 
ELSE  IF(IFLAGl.EQ.l)  THEN 

S0RCE=((UJl^-UJ^)/5.-(UJl*^-UJ>^)*DNl/4.  +  (UJl**3-UJ**3)*DN 
@  2/3.-(UJl**2-UJ**2)*DN3/2.)/DEN0M 
END  IF 
C 

500  RETURN 
END 
C 
C 
C 
C 

C     SUBROUTINE  TO  SOLVE  THE  DISCRETE  ORDINATES  EQUATIONS. 
C 

SUBROUTINE  DO(DEL, EPS, SIGT. CAP, NMAX, ITER) 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  FP(200.50),F1(200,50),F2(200,50).X(50),Z1(50),Z2(50) 
C0MM0N/B1/  F(200,50),G(200.50),Q(200.50),U(50),A(50,50) 
C0MM0N/B2/  M.NR,NT,IS.WIDTH,ERR1,ERR2,ERR3 
MM=M-I 
NTOT=NR+NT 
NR1=NR+1 
C 

XX=DEL*SIGT/2.D0 
SIGT2=SIGT/2.D0 
C 

DO  1  I=1,NT0T 

X(I)=(1.E»-XX/U(I))/(1.D0+XX/U(I)) 
Zl  ( I  )=SIGT2+U(  I  )/DEL 
Z2(I)=SIGT2-U(I)/DEL 

1  CONTINUE 
C 

ITER=0 
100  ITER=ITER+1 
DO  2  I=1,NT0T 
DO  2  K=1,M 
FP(K,I)=F(K,I) 

2  CONTINUE 
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C 

C     SWEEP  FOR  FORWARD  DIRECTIONS 

DO  4  I=NR1,NT0T 

DO  5  K=1,MM 

KK=K+1 

F(KK.  I)=X(I)*F(K,  I)+Q(KK,  I)/Z1(I) 

5  CONTINUE 
4  CONTINUE 

C 

C     SWEEP  FOR  BACKWARD  DIRECTIONS 

DO  6  1=1,  NR 

DO  7  L=1.MM 

K=M-L+1 

KK=K-1 

F(KK,I)=F(K.I)/X(I)+Q(K.I)/Z2(I) 

7  CONTINUE 

6  CONTINUE 
C 

C     CALCULATE  CELL-CENTER  FLUXES  USING  DIAMOND-DIFFERENCE  SCHEME 
DO  8  I=1,NT0T 
DO  9  K=2,M 

G(K,  I)=(F(K-1 , 1)+F(K,  I)  )/2.D0 
9  CONTINUE 

8  CONTINUE 
C 

C     CALCULATE  THE  SCATTERING  SOURCE  TERM 
CAPP=1.D0+CAP 
DO  10  I=1,NT0T 
DO  11  K=2.M 
QS=O.DO 
DO  12  J=1.NT0T 
QS=QS+C(K,J)*A(J,I) 
12  CONTINUE 

Q(K,I)=Q(K,I)+CAPP*(QS-Q(K,I)) 
11  CONTINUE 
10  CONTINUE 
C 

C     CHECK  FOR  CONVERGENCE 
C 

CALL  C0NVRG(F,FP.M,NR,NT.ITER,EPS,ERR1,ERR2,ERR3) 
C 

WRITE(9,910)  ITER,ERR1,ERR2,ERR3 
910  F0RMAT(5X.I2,5X,3(D20.6)) 
C 

C     IF  CONVERGENCE  REQUIREMENTS  ARE  SATISFIED,  RETURN  TO  MAIN 
C     PROGRAM.   IF  NOT,  CONTINUE  WITH  ITERATIONS. 
C 

IF(ERR1.LE. EPS. AND. ERR2.LE.EPS. AND. ERR3.LE.EPS)  GO  TO  200 
C 

C     ADJUST  CONVERGENCE  ACCELERATION  PARAMETER 
150  IF(CAP.GE.2.D-1)  GO  TO  160 

CAP=CAP*1 . 1 
160  CONTINUE 
C 
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C     IF  MAXIMUM  NUMBER  OF  ITERATIONS  HAVE  FAILED  TO  PRODUCE 
C     CONVERGENCE,  CEASE  CALCULATIONS  AND  EXIT  TO  MAIN  PROGRAM. 
C 

IF(ITER.GE.NMAX)  GO  TO  200 
GO  TO  100 
200  RETURN 
END 
C 
C 
C     SUBROUTINE  TO  CHECK  CONVERGENCE  OF  ANGULAR  FLUXES 

SUBROUTINE  CONVRG(F,FP.M,NR,NT,  ITER, EPS, ERR1  ,ERR2,ERR3) 
IMPLICIT  REAL*8  (A-H.O-Z) 
DIMENSION  F(200,50),FP(200,50) 
C 

IF(ITER.EQ.l)  THEN 
ERR1=1. 
ERR2=1. 
ERR3=1. 
GO  TO  500 
END  IF 
C 

C     CHECK  REFLECTED  FLUXES  FOR  CONVERGENCE 
ERR1=0. 
DO  10  1=1. NR 

IF(F(l,I).NE.O.)  ERR1=DMAX1(ERR1,DABS((F(1,I)-FP(1,I))/F(1,I))) 
10  CONTINUE 
C 

C     CHECK  TRANSMITTED  FLUXES  FOR  CONVERGENCE 
ERR2=0. 
NR1=NR+1 
NTOT=NR+NT 
DO  20  I=NR1,NT0T 

IF(F(M,I).NE.O.)  ERR2=DMAX1(ERR2,DABS((F(M,I)-FP(M,I))/F(M,I)}) 
20  CONTINUE 
C 

IF(ERR1. GE.EPS.0R.ERR2.GE.EPS)  THEN 
ERR3=1. 
GO  TO  500 
END  IF 
C 

C     IF  REFLECTED  AND  TRANSMITTED  FLUXES  HAVE  CONVERGED,  CHECK  ALL 
C     FLUXES  FOR  CONVERGENCE 
ERR3=0. 

DO  30  I=1,NT0T 
DO  30  K=1,M 

IF(F(K,I).NE.O.)  ERR3=DMAX1(ERR3,DABS((F(K,I)-FP(K,I))/F(K,I))) 
30  CONTINUE 
C 
C 

500  RETURN 
END 
C 
C 
C     SUBROUTINE  TO  PRINT  OUT  PROBLEM  PARAMETERS  AND  RESULTS 
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SUBROUTINE  OUTPUT(  IFLAG1 ,  IFLAG2 .  ITER , NMAX, EPS ,  SIGT , CC ,  NEXP ,  HINT) 

IMPLICIT  REAL*8(A-H,0-Z) 

REAL*8  JR.JT.JTU.JTOT 

COMMON/Bl/  F(200,50),G(200.50),Q(200,50),U(50),A(50,50) 

C0MM0N/B2/  M,NR,NT,IS,WIDTH,ERR1,ERR2,ERR3 

NT0T=NR+NT 
C 
C 
C     WRITE  OUT  PROBLEM  PARAMETERS 

WRITE(6,900) 

IF(IFLAG2.EQ.0)  THEN 
WRITE(6,901)  WIDTH 

ELSE 

WRITE(6.902)  WIDTH 

END  IF 

WRITE(6,903)  M 

WRITE(6,920)  CC 

WRITE(6,904)  NR 

WRITE(6,905)  NT 

IF(IFLAGl.EQ.O)  THEN 
WRITE(6,906)  U(IS) 

ELSE 

WRITE(6,907)  U(IS) 

END  IF 

WRITE(6,908)  EPS 

WRITE(6,909)  ERR3 

WRITE(6,910)  ITER 
C 
C     CHECK  TO  SEE  IF  CONVERGENCE  WAS  REACHED 

IF(ERR3.LE.EPS)  GO  TO  50 
C     IF  THE  REQUIRED  CONVERGENCE  WAS  NOT  REACHED,  WRITE  OUT 
C     THE  NUMBER  OF  ITERATIONS. 
30  WRITE(6,911)  NMAX 

WRITE(9,911)  NMAX 

GO  TO  100 
C 
C     WRITE  OUT  PROBLEM  RESULTS 

50  CONTINUE 
C 

C     CALL  SUBROUTINE  TO  CALCULATE  THE  REFLECTED  AND  TRANSMITTED 
C     CURRENTS. 

CALL  CURENT (SIGT, JR.JT.JTU.JTOT, NEXP, N INT) 
C 

C     CALCULATE  THE  SCATTERED  FLUX  IN  THE  SOURCE  DIRECTION  BY 
C     SUBTRACTING  THE  UNSCATTERED  FLUX  COMPONENT. 
C 

F(M,IS)=F(M,IS)-F(1,IS)*DEXP(-1.*WIDTH*SIGT/U(IS)) 
C 

WRITE(6,912) 
WRITE(6,913) 
DO  60  I=1,NT0T 

WRITE(6,914)  U(I),F(1,I),M.U(I).F(M,I) 
60  CONTINUE 
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WRITE(6,916)  JR 

WRITE(6,917)  JT 

WRITE(6,918)  JTU 

WRITE(6,919)  JTOT 

IF(NEXP.EQ.l)  THEN 
WRITE(6,921) 

ELSE  IF(NEXP.EQ.2)  THEN 
WRTTE(6,922) 

ELSE  IF(NEXP.EQ.3)  THEN 
WRITE(6.923) 

END  IF 
C 

WRITE(6,915) 
C 
C     OUTPUT  FORMATS 

900  FORMAT( " 1 ' ,2X, 'PROBLEM  PARAMETERS: ' ,/) 

901  FORMAT(5X. 'SLAB  WIDTH      =  \F7.4,  '  CM.') 

902  FORMAT(5X. 'SLAB  WIDTH      =  '  ,F7.4,  '  M.F.P.') 

903  FORMAT(5X. 'NUMBER  OF  NODES  =  ' , 13) 

904  F0RMAT(/,5X,I2, '  REFLECTED  QUADRATURE  DIRECTIONS') 

905  F0RMAT(5X,I2,  '  TRANSMITTED  QUADRATURE  DIRECTIONS') 

906  F0RMAT(5X, 'UNIT  FLUX  INCIDENT  AT  \F9.6,/) 

907  F0RMAT(5X, 'UNIT  CURRENT  INCIDENT  AT  '  .F9.6,/) 

908  F0RMAT(5X,  'CONVERGENCE  CRITERIA  FOR  ANGULAR  FLUXES  =  '  .E11.4) 

909  F0RMAT(5X.  'THE  ANGULAR  FLUX  AT  ALL  POINTS  CONVERGED  TO  WITHIN', Ell 
@.4) 

910  F0RMAT(5X, 'AFTER  ',13,'  ITERATIONS') 

911  F0RMAT(//,5X,'N0  CONVERGENCE  REACHED  AFTER  ',13,'  ITERATIONS') 

912  F0RMAT(4(/),3X,  'REFLECTED  ANGULAR  FLUXES',  17X.  'TRANSMITTED  ANGULAR 
@  FLUXES') 

913  F0RMAT(1X,29('-'),12X.31('-')) 

914  F0RMAT(1X,'F(1,'.F9.6.')  =  ' .E12.5. 12X. 'F( ' , 13, ' , ' ,F9.6, ' )  =  ' ,E12 
@.5) 

915  FORMAT('l') 

916  F0RMAT(////,5X, "THE  REFLECTED  CURRENT  =',F9.6) 

917  F0RMAT(5X,  'THE  COLLIDED  TRANSMITTED  CURRENT   =  '  ,F9.6) 

918  F0RMAT(5X.  'THE  UNCOLLIDED  TRANSMITTED  CURRENT  =  \F9.6) 

919  F0RMAT(/,5X,'THE  TOTAL  ESCAPING  CURRENT        =  \F9.6) 

920  F0RMAT(5X, 'C  =  '  ,F7.4) 

921  F0RMAT(3(/),2X,  'APPROXIMATED  SCATTERING  KERNEL  METHOD  USED  (LINEAR 
@  FLUX  EXPANSION)  '  ) 

922  F0RMAT(3(/),2X,  'APPROXIMATED  SCATTERING  KERNEL  METHOD  USED  (QUADRA 
OTIC  FLUX  EXPANSION)') 

923  F0RMAT(3(/),2X,  'APPROXIMATED  SCATTERING  KERNEL  METHOD  USED  (CUBIC 
@FLUX  EXPANSION)') 

100  CONTINUE 

RETURN 

END 
C 
C 
C 
C     SUBROUTINE  TO  CALCULATE  ESCAPING  CURRENTS 

SUBROUTINE  CURENT(SICT, JR. JT, JTU,  JTOT, NEXP.N INT) 

IMPLICIT  REAL*8  (A-H.O-Z) 
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REAL*8  JR.JT.JTT.JTU.JTOT.NORM 

DIMENSION  FF(50),UU(50) 

COMMON/B1/  F(200,50),G(200.50),Q(200,50),U(50),A(50,50) 

C0MM0N/B2/  M,NR,NT.IS,WIDTH,ERR1,ERR2.ERR3 
C 

NR1=NR+1 

NTOT=NR+NT 
C 
C     CALCULATE  REFLECTED  CURRENT 

JR=0. 

DO  100  1=1,  NR 

CALL  N0RMAL(U,  1,1, NORM, NEXP.NINT) 

JR=JR+F(1.I)*N0RM 
100  CONTINUE 
C 

JK=-1.*JR 
C     CALCULATE  TOTAL  TRANSMITTED  CURRENT 

JTT=0. 

DO  200  I=NR1,NT0T 

CALL  NORMAL(U,  1,1, NORM, NEXP.NINT) 

JTT=JTT+F(M.  I)*NORM 
200  CONTINUE 
C 
C     CALCULATE  THE  COLLIDED  TRANSMITTED  CURRENT 

DO  210  I=NR1,NT0T 

FF(I)=F(M,I) 
210  CONTINUE 

FF( IS)=F(M,  IS)-F(  1 ,  IS)*DEXP(-1 . *WIDTH*SIGT/U( IS)  ) 
C 

JT=0. 

DO  220  I=NR1.NT0T 

CALL  NORMAL(U, 1,1, NORM, NEXP.NINT) 

JT=JT+FF(I)*NORM 
220  CONTINUE 
C 

C  CALCULATE  THE  UNCOLLIDED  TRANSMITTED  CURRENT  BY  SUBTRACTING  THE 
C  COLLIDED  TRANSMITTED  CURRENT  FROM  THE  TOTAL  TRANSMITTED  CURRENT 
C 

JTU=JTT-JT 
C 

JTOT=JR+JTT 
C 

RETURN 

END 
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ABSTRACT 

A  fundamental  difficulty  with  discrete-ordinates  transport 
calculations  is  the  accurate  evaluation  of  the  scattering  source  term 
for  highly  anisotropic  problems  (i.e.,  those  problems  that  are 
characterized  by  anisotropic  fluxes  as  well  as  by  anisotropic 
scattering).  For  such  problems,  the  conventional  Legendre  expansion  of 
the  differential  scattering  cross  section  often  yields  very  inaccurate 
and  even  negative  angular  fluxes.  Alternatively,  the  use  of  the  exact 
cross  sections  for  transfer  from  one  discrete  direction  to  another  (the 
"exact  kernel  method")  provides  better  accuracy,  but  suffers  from  an 
angular  redistribution  problem  when  low-order  quadrature  is  used. 

In  this  work,  a  new  method  for  evaluating  the  scattering  source 
term  is  proposed.  This  new  method  represents  both  the  angular  flux  and 
the  scattering  cross  section  by  piecewise  polynomial  expansions  so  that 
the  resulting  scattering  source  term  can  be  integrated  analytically. 
The  integrated  results  can  be  used  in  a  standard  discrete-ordinates  code 
with  minor  modification  of  the  scattering  source  term  calculation. 

From  the  results  of  several  slab  albedo  transport  problems,  it  is 
apparent  that  the  new  method  provides  better  accuracy  than  does  the 
Legendre  expansion  method,  and  also  eliminates  the  angular 
redistribution  problem  which  limits  the  exact  kernel  method  for  low 
quadrature  orders.  The  accuracy  which  can  be  obtained  with  the  new 
method,  however,  depends  on  the  degree  of  anisotropy  in  the  angular 
flux,  the  number  of  discrete  ordinates  used,  and  the  order  of  the 
piecewise  polynomial  expansion  used  to  represent  the  angular  flux. 
Often  the  simple  linear  approximation  gives  the  best  results  for  highly 
anisotropic  problems,  since  this  approximation  is  the  only  one 
guaranteed  to  always  yield  a  positive  source  term. 


